Merge pull request #33 from tjjarvinen/isolated_system_support

Implement isolated system for AB (#32)
JuliaMolSim · Jan 13, 2025 · 1bdd344 · 1bdd344
2 parents 71fd91d + 9d3365d
commit 1bdd344
Show file tree

Hide file tree

Showing 3 changed files with 52 additions and 3 deletions.
diff --git a/src/atoms_base.jl b/src/atoms_base.jl
@@ -4,13 +4,34 @@ using Unitful
 
 
 function PairList(ab::AtomsBase.AbstractSystem, cutoff::Unitful.Length; length_unit=unit(cutoff))
-    cell = ustrip.(length_unit, hcat( cell_vectors(ab)... )' )
-    pbc = periodicity(ab)
     r = map( 1:length(ab) ) do i
         # Need to have SVector here for PairList to work
         # if position does not give SVector
         SVector( ustrip.(length_unit, position(ab,i))...)
     end
-    nlist = PairList(r, ustrip(length_unit, cutoff), cell, pbc; int_type=Int)
+    tmp_cell = cell(ab)
+    if tmp_cell isa AtomsBase.IsolatedCell{3, <:Any}
+        # this only works in 3D
+        max_x = maximum(q->q[1], r)
+        max_y = maximum(q->q[2], r)
+        max_z = maximum(q->q[3], r)
+        min_x = minimum(q->q[1], r)
+        min_y = minimum(q->q[2], r)
+        min_z = minimum(q->q[3], r)
+        cell_matrix = SMatrix{3,3, typeof(max_x), 9}(
+            # add +1 for cell size
+            max_x - min_x + 1, 0., 0.,
+            0., max_y - min_y + 1, 0.,
+            0., 0., max_z - min_z + 1,
+        )
+    elseif tmp_cell isa AtomsBase.IsolatedCell
+        D = n_dimensions(ab)
+        throw( error("NeighbourLists.jl does not support $D-dimensional isolated AtomsBase systems yet.") )
+    else
+        cell_matrix = ustrip.(length_unit, hcat( cell_vectors(ab)... )' )
+    end
+    pbc = periodicity(ab)
+
+    nlist = PairList(r, ustrip(length_unit, cutoff), cell_matrix, pbc; int_type=Int)
     return nlist
 end
diff --git a/test/Project.toml b/test/Project.toml
@@ -1,4 +1,5 @@
 [deps]
+AtomsBase = "a963bdd2-2df7-4f54-a1ee-49d51e6be12a"
 AtomsBuilder = "f5cc8831-eeb7-4288-8d9f-d6c1ddb77004"
 Distances = "b4f34e82-e78d-54a5-968a-f98e89d6e8f7"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"

diff --git a/test/test_atoms_base.jl b/test/test_atoms_base.jl
@@ -1,4 +1,5 @@
 #using ASEconvert
+using AtomsBase
 using AtomsBuilder
 using NeighbourLists
 using Test
@@ -12,6 +13,32 @@ using Unitful
     id, r = neigs(nlist,1)
 
     @test length(id) == 12
+
+    @testset "AtomsBase Isolated System" begin
+        isys = isolated_system([
+            :H => [0., 0., 0.]u"Å",
+            :H => [0., 0., 2.]u"Å",
+            :H => [0., 10., 0.]u"Å"
+        ])
+        nlist = PairList(isys, 5.0u"Å")
+        id, r = neigs(nlist,1)
+        @test length(id) == 1
+        @test id[1] == 2
+        @test all( r[1] .≈ [0., 0., 2.] )
+        id, r = neigs(nlist,2)
+        @test length(id) == 1
+        @test id[1] == 1
+        @test all( r[1] .≈ [0., 0., -2.] )
+        id, r = neigs(nlist,3)
+        @test length(id) == 0
+
+        isys2d = isolated_system([
+            :H => [ 0., 0.]u"Å",
+            :H => [ 0., 2.]u"Å",
+            :H => [ 10., 0.]u"Å"
+        ])
+        @test_throws ErrorException PairList(isys2d, 5.0u"Å")
+    end
 end