@testset "SOAP Descriptor" begin

@info("Testing SOAP Descriptor for DC Si with cutoff=6.5, n=4, l=4")
using JuLIP
using QUIP, Test

at = bulk(:Si, cubic=true)
desc = soap(at, 6.5, n_max=4, l_max=4)
# Reference is from the output of the QUIP command as follows
# quip atoms_filename=test.xyz descriptor_str="soap cutoff=6.5 l_max=4 n_max=4 atom_sigma=0.5 normalise=T n_Z=1 Z={14} n_species=1 species_Z={14}"
soap_ref = [ 0.23864044377742688E+000,  0.64340429976653843E-038,  0.92117979653018819E-038, 
             0.48626301215687845E-005,  0.18876337954462125E-005,  0.29409765463423732E+000, 
            -0.11518886345720650E-036,  0.11050218320476954E-035,  0.11763958794274931E-002,  
             0.45589948982910106E-003,  0.18122123201805027E+000,  0.11109715229581757E-033,  
             0.20740200612211439E-033,  0.14230028097299613E+000,  0.55054202072415959E-001,  
             0.51605689877931338E+000, -0.13873220359922910E-035,  0.53089215594650063E-036, 
            -0.50839942118860665E-003, -0.50804299476442307E-003,  0.44970747195745586E+000,  
             0.56664079661428650E-034, -0.62628772367525336E-034, -0.86970568361713990E-001, 
            -0.86763453783360533E-001,  0.55798321223812897E+000,  0.17093065575762176E-033,  
             0.32536812908176952E-031,  0.26577177885526940E-001,  0.68368050294466742E-001,  
             0.37916397221481701E-002,  0.29040731747973391E-037, -0.11050630810221818E-037,  
             0.11135545482166505E-004,  0.46298939173893863E-005,  0.33041486666568242E-002, 
            -0.88957883103660356E-036, -0.32224333298567884E-036,  0.19049288398824812E-002,  
             0.79069210886284307E-003,  0.57978317308968498E-002, -0.67817971363747986E-035,  
             0.24028731016100991E-034, -0.82324728991149839E-003, -0.88112752017744061E-003,  
             0.30121742054713143E-004,  0.23815358022366506E-036,  0.18522830192896921E-037,  
             0.12750339845444976E-004,  0.56779863069817660E-005]
soap_now = vcat(desc[1,:]...)
println("QUIP.jl SOAP:",soap_now)
println("Reference:",soap_ref)
println(@test soap_now  ≈  soap_ref)

end
