ACE Descriptors

This tutorial demonstrates a simple use of ACE descriptors.

using ACEpotentials, MultivariateStats, Plots

Load a tiny silicon dataset, which has the isolated atom, 25 diamond-like (dia) configurations, 25 beta-tin-like (bt) configurations, and 2 liquid (liq) configurations.

dataset, _, _ = ACEpotentials.example_dataset("Si_tiny");

An ACE basis specifies a vector of invariant features of atomic environments and can therefore be used as a general descriptor.

basis = ACE1x.ace_basis(elements = [:Si],
                        rcut = 5.5,
                        order = 3,        # body-order - 1
                        totaldegree = 8);

Compute an averaged structural descriptor for each configuration.

descriptors = []
config_types = []
for atoms in dataset
    struct_descriptor = sum(site_descriptors(basis, atoms)) / length(atoms)
    push!(descriptors, struct_descriptor)
    push!(config_types, atoms.data["config_type"].data)
end

Finally, extract the descriptor principal components and plot. Note the segregation by configuration type.

descriptors = hcat(descriptors...)
M = fit(PCA, descriptors; maxoutdim=3, pratio=1)
descriptors_trans = transform(M, descriptors)
p = scatter(
    descriptors_trans[1,:], descriptors_trans[2,:], descriptors_trans[3,:],
    marker=:circle, linewidth=0, group=config_types, legend=:right)
plot!(p, xlabel="PC1", ylabel="PC2", zlabel="PC3", camera=(20,10))
Example block output

The basis used above uses defaults that are suitable for regression of a potential energy surface, but other defaults might be better when using the ACE descriptor for other tasks such as classification. The following short script shows how to make some changes of this kind:

model = acemodel(elements = [:Si,], order = 3, totaldegree = 10,
       pair_transform = (:agnesi, 1, 4, 0.0),
       pair_envelope = (:x, 0, 2),
       transform = (:agnesi, 1, 4, 0.0),
       envelope = (:x, 0, 2),
       r0 = :bondlen, # default, could specify explicitly
       )
basis = model.basis
JuLIP.MLIPs.IPSuperBasis{JuLIP.MLIPs.IPBasis}(JuLIP.MLIPs.IPBasis[PolyPairBasis{ACE1.OrthPolys.TransformedPolys{Float64, ACE1.Transforms.Normalized{Float64, ACE1.Transforms.Agnesi2Transform{Float64, Int64}}, ACE1.OrthPolys.OrthPolyBasis{Float64}, ACE1.OrthPolys.OneEnvelope}, 1}(JuLIP.Potentials.SZList{1}(AtomicNumber[<14>]), ACE1.OrthPolys.TransformedPolys{Float64, ACE1.Transforms.Normalized{Float64, ACE1.Transforms.Agnesi2Transform{Float64, Int64}}, ACE1.OrthPolys.OrthPolyBasis{Float64}, ACE1.OrthPolys.OneEnvelope}[ACE1.OrthPolys.TransformedPolys{Float64, ACE1.Transforms.Normalized{Float64, ACE1.Transforms.Agnesi2Transform{Float64, Int64}}, ACE1.OrthPolys.OrthPolyBasis{Float64}, ACE1.OrthPolys.OneEnvelope}(ACE1.OrthPolys.OrthPolyBasis{Float64}(0, 0.0, 2, 1.0, [2.2360689091951995, 7.099309346735727, 5.291522024409, 4.738063448424564, 4.484386643542148, 4.344482007955068, 4.258391894831597, 4.201377361610357, 4.1615577274183035, 4.132600058800379], [0.0, -1.1832206894942505, -1.7638467747583702, -1.8952356846820837, -1.943249621114157, -1.965382380919459, -1.9771390388551457, -1.9840203025885013, -1.9883453766477837, -1.9912174403384049], [0.0, 0.0, -0.7453572968815695, -0.8954065440091865, -0.9464598126125228, -0.9688018347417563, -0.9801840327648208, -0.9866112526443501, -0.9905222428825616, -0.9930416275551956], [0.0005, 0.0015, 0.0025, 0.0035, 0.0045, 0.0055, 0.0065, 0.0075, 0.0085, 0.0095  …  0.9905, 0.9915, 0.9925, 0.9935, 0.9945, 0.9955, 0.9965, 0.9975, 0.9985, 0.9995], [0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001  …  0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001]), ACE1.Transforms.Normalized{Float64, ACE1.Transforms.Agnesi2Transform{Float64, Int64}}(ACE1.Transforms.Agnesi2Transform{Float64, Int64}(2.4, 1, 4, 0.2727272727272727, 0.0), 1.0, 0.6094563632576547, 6.0), 0.0, 6.0, ACE1.OrthPolys.OneEnvelope());;], [0;;]), RPIBasis{Float64, BasicPSH1pBasis{Float64, 1, ACE1.OrthPolys.TransformedPolys{Float64, ACE1.Transforms.MultiTransform{1, ACE1.Transforms.AffineT{Float64, ACE1.Transforms.Normalized{Float64, ACE1.Transforms.Agnesi2Transform{Float64, Int64}}}}, ACE1.OrthPolys.OrthPolyBasis{Float64}, ACE1.OrthPolys.OneEnvelope}}, 1, ACE1.DAGEvaluator}(PIBasis{BasicPSH1pBasis{Float64, 1, ACE1.OrthPolys.TransformedPolys{Float64, ACE1.Transforms.MultiTransform{1, ACE1.Transforms.AffineT{Float64, ACE1.Transforms.Normalized{Float64, ACE1.Transforms.Agnesi2Transform{Float64, Int64}}}}, ACE1.OrthPolys.OrthPolyBasis{Float64}, ACE1.OrthPolys.OneEnvelope}}, 1, ACE1.DAGEvaluator}(BasicPSH1pBasis{Float64, 1, ACE1.OrthPolys.TransformedPolys{Float64, ACE1.Transforms.MultiTransform{1, ACE1.Transforms.AffineT{Float64, ACE1.Transforms.Normalized{Float64, ACE1.Transforms.Agnesi2Transform{Float64, Int64}}}}, ACE1.OrthPolys.OrthPolyBasis{Float64}, ACE1.OrthPolys.OneEnvelope}}(ACE1.OrthPolys.TransformedPolys{Float64, ACE1.Transforms.MultiTransform{1, ACE1.Transforms.AffineT{Float64, ACE1.Transforms.Normalized{Float64, ACE1.Transforms.Agnesi2Transform{Float64, Int64}}}}, ACE1.OrthPolys.OrthPolyBasis{Float64}, ACE1.OrthPolys.OneEnvelope}(ACE1.OrthPolys.OrthPolyBasis{Float64}(0, -1.0, 2, 1.0, [0.5590172272987999, 3.5496546733678636, 2.6457610122045, 2.3690317242122805, 2.2421933217710746, 2.1722410039775335, 2.129195947415799, 2.1006886808051783, 2.0807788637091518, 2.066300029400189, 2.055429679150542, 2.047054944278475, 2.040463473582576], [0.0, 2.3664339838736135, 0.881914237446129, 0.4737960395301985, 0.29894370065691694, 0.20685862305807545, 0.15205690856065218, 0.1166683782166776, 0.09243348706136743, 0.07508258906178471, 0.0622191370400381, 0.05240996282372767, 0.044753712522501335, 0.03865970882897952], [0.0, 0.0, -0.7453572968815678, -0.8954065440091885, -0.9464598126125205, -0.9688018347417586, -0.9801840327648191, -0.9866112526443519, -0.9905222428825597, -0.9930416275551971, -0.9947392198156224, -0.9959255551493605, -0.9967800225810635, -0.997411318554425], [-0.999, -0.997, -0.995, -0.993, -0.991, -0.989, -0.987, -0.985, -0.983, -0.981  …  0.981, 0.983, 0.985, 0.987, 0.989, 0.991, 0.993, 0.995, 0.997, 0.999], [0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001  …  0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001]), ACE1.Transforms.MultiTransform{1, ACE1.Transforms.AffineT{Float64, ACE1.Transforms.Normalized{Float64, ACE1.Transforms.Agnesi2Transform{Float64, Int64}}}}(JuLIP.Potentials.SZList{1}(AtomicNumber[<14>]), ACE1.Transforms.AffineT{Float64, ACE1.Transforms.Normalized{Float64, ACE1.Transforms.Agnesi2Transform{Float64, Int64}}}[ACE1.Transforms.AffineT{Float64, ACE1.Transforms.Normalized{Float64, ACE1.Transforms.Agnesi2Transform{Float64, Int64}}}(ACE1.Transforms.Normalized{Float64, ACE1.Transforms.Agnesi2Transform{Float64, Int64}}(ACE1.Transforms.Agnesi2Transform{Float64, Int64}(2.4, 1, 4, 0.2727272727272727, 0.0), 1.0, 0.6094563632576547, 6.0), 0.0, 1.0, -1.0, 1.0);;]), 0.0, 6.0, ACE1.OrthPolys.OneEnvelope()), ACE1.SphericalHarmonics.SHBasis{Float64}(2, ACE1.SphericalHarmonics.ALPCoefficients{Float64}([6.9280545871776e-310, 6.92819467604806e-310, 0.0, 1.9364916731037085, 0.0, 0.0], [6.9281911422036e-310, 6.9281911422036e-310, 6.9281911422036e-310, -0.5773502691896257, 6.9281911422036e-310, 6.9281911422036e-310])), JuLIP.Potentials.SZList{1}(AtomicNumber[<14>]), ACE1.RPI.PSH1pBasisFcn[znlm[0|1,0,0], znlm[0|1,1,-1], znlm[0|1,1,0], znlm[0|1,1,1], znlm[0|1,2,-2], znlm[0|1,2,-1], znlm[0|1,2,0], znlm[0|1,2,1], znlm[0|1,2,2], znlm[0|2,0,0]  …  znlm[0|6,0,0], znlm[0|6,1,-1], znlm[0|6,1,0], znlm[0|6,1,1], znlm[0|7,0,0], znlm[0|8,0,0], znlm[0|9,0,0], znlm[0|10,0,0], znlm[0|11,0,0], znlm[0|12,0,0]], UnitRange{Int64}[1:45;;]), JuLIP.Potentials.SZList{1}(AtomicNumber[<14>]), (ACE1.InnerPIBasis([1, 1, 1, 1, 1, 1, 1, 1, 1, 1  …  3, 3, 3, 3, 3, 3, 3, 3, 3, 3], [1 0 0; 10 0 0; … ; 19 19 19; 19 19 28], Dict{ACE1.PIBasisFcn, Int64}(ACE1.PIBasisFcn{3, ACE1.RPI.PSH1pBasisFcn}(<14>, (znlm[14|1,1,0], znlm[14|1,2,-1], znlm[14|2,1,1]), ACE1.RPI.PSH1pBasisFcn) => 166, ACE1.PIBasisFcn{2, ACE1.RPI.PSH1pBasisFcn}(<14>, (znlm[14|1,0,0], znlm[14|9,0,0]), ACE1.RPI.PSH1pBasisFcn) => 21, ACE1.PIBasisFcn{3, ACE1.RPI.PSH1pBasisFcn}(<14>, (znlm[14|1,1,0], znlm[14|2,1,0], znlm[14|3,0,0]), ACE1.RPI.PSH1pBasisFcn) => 172, ACE1.PIBasisFcn{3, ACE1.RPI.PSH1pBasisFcn}(<14>, (znlm[14|1,0,0], znlm[14|1,1,0], znlm[14|2,1,0]), ACE1.RPI.PSH1pBasisFcn) => 101, ACE1.PIBasisFcn{2, ACE1.RPI.PSH1pBasisFcn}(<14>, (znlm[14|1,0,0], znlm[14|8,0,0]), ACE1.RPI.PSH1pBasisFcn) => 20, ACE1.PIBasisFcn{3, ACE1.RPI.PSH1pBasisFcn}(<14>, (znlm[14|1,0,0], znlm[14|1,2,-1], znlm[14|2,2,1]), ACE1.RPI.PSH1pBasisFcn) => 112, ACE1.PIBasisFcn{3, ACE1.RPI.PSH1pBasisFcn}(<14>, (znlm[14|1,1,0], znlm[14|2,0,0], znlm[14|2,1,0]), ACE1.RPI.PSH1pBasisFcn) => 169, ACE1.PIBasisFcn{2, ACE1.RPI.PSH1pBasisFcn}(<14>, (znlm[14|1,0,0], znlm[14|4,0,0]), ACE1.RPI.PSH1pBasisFcn) => 16, ACE1.PIBasisFcn{2, ACE1.RPI.PSH1pBasisFcn}(<14>, (znlm[14|2,1,0], znlm[14|2,1,0]), ACE1.RPI.PSH1pBasisFcn) => 63, ACE1.PIBasisFcn{2, ACE1.RPI.PSH1pBasisFcn}(<14>, (znlm[14|2,1,-1], znlm[14|4,1,1]), ACE1.RPI.PSH1pBasisFcn) => 61…), Dict{Any, Int64}(), 1:206, <14>, ACE1.DAG.CorrEvalGraph{Int64, Int64}([(1, 0), (2, 0), (3, 0), (4, 0), (5, 0), (6, 0), (7, 0), (8, 0), (9, 0), (10, 0)  …  (10, 79), (10, 81), (73, 19), (73, 28), (73, 32), (74, 28), (75, 19), (78, 19), (83, 19), (83, 28)], [1, 0, 0, 0, 0, 0, 0, 0, 0, 2  …  197, 198, 199, 200, 201, 202, 203, 204, 205, 206], 45, 95)),), ACE1.DAGEvaluator()), (sparse([1, 2, 3, 10, 11, 12, 16, 17, 18, 19  …  93, 93, 94, 95, 96, 97, 98, 98, 99, 100], [1, 1, 1, 1, 1, 1, 1, 1, 1, 1  …  197, 198, 199, 200, 201, 202, 203, 204, 205, 206], [-0.6371868657466305, 0.13327692844605593, -0.016791269897827976, -1.1036400253887322, 0.2308424115452905, -0.029083332586639907, 1.424793146179496, -0.2980162718377561, 0.03754642092008926, -0.44602954226730596  …  0.5773502691896257, -0.5773502691896257, 1.4142135623730951, 1.0, 1.0, 1.4142135623730951, -1.6329931618554518, 0.8164965809277259, 2.449489742783178, 1.4142135623730951], 100, 206),), (1:100,))])
  • [pair_]transfor = (:agnesi, 1, 4, 0.0) : this generates a transform that behaves as t' ~ r^3 near zero and t' ~ r^-2 near the cutoff, but then actually enforces t' = 0 at the cutoff.
  • [pair_]envelope = (:x, 0, 2) : this generates an envelope that is ~ (x - xcut)^2 at the cutoff and just ~ 1 for r -> 0.

This page was generated using Literate.jl.