Merge pull request #52 from cesmix-mit/sw/change_descr_output

Update ACE force descriptor output

Project.toml

@@ -1,7 +1,7 @@
 name = "InteratomicPotentials"
 uuid = "a9efe35a-c65d-452d-b8a8-82646cd5cb04"
 authors = ["Dallas Foster <[email protected]>"]
-version = "0.2.7"
+version = "0.2.8"
 
 [deps]
 ACE1 = "e3f9bc04-086e-409a-ba78-e9769fe067bb"

README.md

@@ -12,6 +12,8 @@ This package is part of the CESMIX molecular modeling suite. This package is als
 
 This package is a work in progress. 
 
+**WARNING**: As of v0.2.8, SNAP implementation is inconsistent with ACE implementation, and additionally, may produce incorrect values! (See [this issue](https://github.com/cesmix-mit/AtomisticComposableWorkflows/issues/10)).  Use at your own risk!
+
 ## Working Example
 
 In order to compute the interatomic energy of a system, or the forces between atoms in a system, the user has to

src/BasisSystems/ACE/ace.jl

@@ -77,11 +77,12 @@ end
 
 function compute_force_descriptors(A::AbstractSystem, ace::ACE)
     ftemp = ACE1.forces(ace.rpib, convert_system_to_atoms(A))
-    f = [zeros(3, length(ace)) for i = 1:length(A)]
+    # TODO: settle on efficient and consistent descriptor layout
+    f = [ [zeros(length(ace)) for _ in 1:3] for i in 1:length(A)]
     for i = 1:length(A)
         for j = 1:3 
             for k = 1:length(ace)
-                f[i][j, k] = ftemp[k][i][j]
+                f[i][j][k] = ftemp[k][i][j]
             end
         end
     end

src/types/empirical_potential.jl

@@ -65,3 +65,13 @@ function virial_stress(s::AbstractSystem, p::EmpiricalPotential)
     end
     SVector{6}(v) * ENERGY_UNIT
 end
+
+function force(s::AbstractSystem, p::EmpiricalPotential) 
+    eandf = energy_and_force(s,p)
+    eandf.f
+end
+
+function potential_energy(s::AbstractSystem, p::EmpiricalPotential) 
+    eandf = energy_and_force(s,p)
+    eandf.e
+end

test/ACE/ace_test.jl

@@ -23,17 +23,28 @@ ace = ACE( species = [:Ar],
            polynomial_degree = 8,
            rcutoff = 2.0)
 @test isa(ace, BasisSystem)
-e = sum(compute_local_descriptors(system, ace))
-@test isa(e, AbstractVector)
-f = compute_local_descriptors(system, ace)
-@test all(isa.(f, (AbstractVector,)))
+e_descrs_l = compute_local_descriptors(system, ace)
+@test all(isa.(e_descrs_l, Vector{<:Real}))
+@test length(e_descrs_l) == 3
+@test length(e_descrs_l[1]) == 8
+e_descrs_g = sum(e_descrs_l)
+@test isa(e_descrs_g, AbstractVector)
+
+f_descrs = compute_force_descriptors(system, ace)
+@test typeof(f_descrs) <: Vector{<:Vector{<:Vector{<:Real}}}
+@test length(f_descrs) == 3 # num of atoms
+@test length(f_descrs[1]) == 3 # x,y,z
+@test length(f_descrs[1][1]) == 8 # number of descriptors
+
 v = compute_virial_descriptors(system, ace)
 @test isa(v, AbstractArray)
 
-print(ace)
+#print(ace)
 lbp = LBasisPotential(ace)
 @test isa(lbp, BasisPotential)
 
 basis = get_rpi(ace);
 
-
+@testset "Reference ACE calculations" begin
+    include("reference_ace_test.jl")
+end

test/ACE/reference_ace_test.jl

@@ -0,0 +1,50 @@
+
+### set up ACE potential
+ace = ACE(species           = [:Ti, :Al],
+          body_order        = 3,
+          polynomial_degree = 6,
+          wL                = 1.5,
+          csp               = 1.0,
+          r0                = 2.9,
+          rcutoff           = 5.5 )
+
+coeffs = [0.025591381519026762, 0.03527125788791906, 0.030612348271342734, 0.06646319273427727, -0.007326117584533097, 0.0021476223879284516, 0.007005564677788114, 0.007769244005502122, 0.0033019143419533176, -0.024266397752728517, -0.03511058549188825, -0.004486106026460792, 0.07649832144216986, -0.017295005584346018, 0.0348519410800185, -0.0026935081045876344, -0.0036996351104090115, 0.04250332320742131, -0.06611126598243479, 0.07452744399669442, -0.08807382022058645, 0.006553101218837121, -0.02825330387435087, -0.005070437887508557, 0.017488241826946662, -0.041461388491636234, -0.050152804966179194, 0.014554551186620662, 0.005494466857846328, 0.03395869840669037, -0.12004390275966798, 0.07758118243125994, 0.024624168020804672, 0.0006581992277555695, -0.002196641935532242, 0.03231551745953874, 0.0005431753297032715, 0.009602374511533056, 0.028907266845791348, 0.03557855646347803, 0.000832998634326787, 0.019238505326450918, -0.007863457928993406, 0.03497657242548427, -0.058485491203844206, -0.025527625067137013, -0.003851837725125408, 0.019472633328804008, -0.04975455754968226, 0.008243807089528446, 0.020612783411412677, -0.07411984524326856, 0.007005564677788615, 0.0077692440055020795, 0.003301914341951156, -0.024266397752728874, -0.03511058549188732, -0.004486106026461139, 0.004050167320520888, 0.011275083723136878, -0.009533633282696359, -0.016652089366136488, 0.005947187981081792, -0.0086386178798077, 0.027556838876613178, -0.01794755394550558, 0.0328518497817209, -0.0444008944069233, -0.04142464521909121, 0.014939466653767677, -0.0013061815492572404, -0.008399904141687925, -0.013070571180237286, 0.07022679858374972, 0.03655463426663164, -0.02425878114877371, -0.013089322405632224, 0.007663504514768707, -0.0006932536563853398, -0.015392489165582057, -0.005333834581033578, 0.000966860983206308, -0.06259246382383571, 0.04896321372972445, -0.012976299346956766, 0.00575471543255263, -0.010710826925328487, 0.009130893987440367, 0.025455356200891677, 0.03737467186835743, -0.061410072131176816, 0.05891873535070835, 0.02179281899408886, 0.02640823532251172, 0.0002904232787473565, -0.028649695323579742, -0.018788426151163752, -0.004911376520223526, 0.034242726688142995, -0.008960425717451335, -0.04627434272845332, 0.05321984323617818, 0.013802856612787245, 0.023560957961937884]
+
+lbp = LBasisPotential(coeffs,[0.0],ace)
+
+
+# set up small nonorthogonal TiAl system 
+atom1 = Atom(:Ti, (@SVector [1.47692848, 1.36039173,1.98890986]) * u"Å")
+atom2 = Atom(:Al, (@SVector [-0.07703613, 0.04783871, -0.03196409]) * u"Å")
+atoms = [atom1, atom2]
+box = [[2.8186357315284174, 0.0, 0.0], 
+       [0.0, 2.8186357315284174, 0.0], 
+       [0.0, 0.1959273402867757, 4.064092262967999]]u"Å"
+bcs = [Periodic(), Periodic(), Periodic()]
+small_nonorth_sys = FlexibleSystem(atoms, box, bcs)
+
+#= 
+TODO: document origin of reference values better 
+    ACE1 v0.12.0
+    JuLIP v0.14.5
+    Julia 1.9.2+0.aarch64.apple.darwin14
+
+    cross-referenced w/ ML-PACE, exported via ACE1pack
+=#
+
+pe = potential_energy(small_nonorth_sys,lbp)
+@test pe ≈ -2.5861943240822214
+
+ref_f1 = [-0.4461309280449358, 0.5403359378738349, -0.1331972765358]
+ref_f2 = [0.4461309280449358, -0.5403359378738349, 0.13319727653580002]
+forces = force(small_nonorth_sys,lbp)
+@test forces[1] ≈ ref_f1
+@test forces[2] ≈ ref_f2 
+(e,f) = energy_and_force(small_nonorth_sys,lbp) 
+@test e == pe
+@test f == forces
+
+
+#TODO: test virial_stress, virial
+#TODO: test functions dispatched off of descriptors 
+#TODO: test get_rcut, get_species, get_params