Infinite energy

src/moves.jl

@@ -1,5 +1,16 @@
 using ComponentArrays
 
+function Arianna.perform_action!(system::Particles, action::Action)
+    e₁, e₂ = perform_action!(system, action)
+    if isinf(e₂)
+        action.δe = zero(typeof(system.energy[1]))
+    else
+        action.δe = e₂ - e₁
+        system.energy[1] += action.δe
+    end
+    return e₁, e₂
+end
+
 ###############################################################################
 # SIMPLE DISPLACEMENT
 
@@ -22,7 +33,7 @@ function update_position!(system::Particles, action::Displacement)
     @inbounds system.position[action.i] = system.position[action.i] + action.δ
 end
 
-function Arianna.perform_action!(system::Particles, action::Displacement)
+function perform_action!(system::Particles, action::Displacement)
     neighbour_list = get_neighbour_list(system)
     e₁ = destroy_particle!(system, action.i, neighbour_list)
     update_position!(system, action)
@@ -32,7 +43,6 @@ function Arianna.perform_action!(system::Particles, action::Displacement)
     end
     e₂ = create_particle!(system, action.i, neighbour_list)
     action.δe = e₂ - e₁
-    system.energy[1] += action.δe
     return e₁, e₂
 end
 
@@ -92,12 +102,11 @@ function update_species_list!(species_list, swap_species, i, j)
     species_list.sp_heads[i], species_list.sp_heads[j] = species_list.sp_heads[j], species_list.sp_heads[i]
 end
 
-function Arianna.perform_action!(system::Particles, action::DiscreteSwap)
+function perform_action!(system::Particles, action::DiscreteSwap)
     i, j = action.i, action.j
     spi, spj = get_species(system, i), get_species(system, j)
     e₁, e₂ = swap_particle_species!(system, spi, i, spj, j)
     action.δe = e₂ - e₁
-    system.energy[1] += action.δe
     update_species_list!(system.species_list, action.species, i, j)
     return e₁, e₂
 end
@@ -157,12 +166,11 @@ mutable struct MoleculeFlip{F<:AbstractFloat} <: Action
     δe::F
 end
 
-function Arianna.perform_action!(system::Particles, action::MoleculeFlip)
+function perform_action!(system::Particles, action::MoleculeFlip)
     i, j = action.i, action.j
     spi, spj = system.species[i], system.species[j]
     e₁, e₂ = swap_particle_species!(system, spi, i, spj, j)
     action.δe = e₂ - e₁
-    system.energy[1] += action.δe
     return e₁, e₂
 end
 

test/molecules_test.jl

@@ -14,25 +14,6 @@ Length = 3
 d = 3
 temperature = 2.0
 density = 1.2
-box = @SVector fill(typeof(temperature)((N / density)^(1 / d)), d)
-position_1 = [box .* @SVector rand(rng, d) for i in 1:Int(N // Length), m in 1:M]
-position = []
-for r in position_1
-    push!(position, r)
-    push!(position, r .+ @SVector [0.1, 0.1, 0.1])
-    push!(position, r .- @SVector [0.1, 0.1, 0.1])
-end
-position = Vector{SVector{3,Float64}}(position)
-molecule = Vector{Int}()
-species = Vector{Int}()
-for i in 1:Int(N // Length)
-    push!(species, 1)
-    push!(species, 2)
-    push!(species, 3)
-    push!(molecule, i)
-    push!(molecule, i)
-    push!(molecule, i)
-end
 
 function create_bond_matrix(N::Int)
     # Create a vector to store the SVectors, each containing a pair of integers
@@ -46,16 +27,15 @@ function create_bond_matrix(N::Int)
     return matrix
 end
 
+chains = load_chains("test/molecule.xyz", args=Dict("temperature" => [temperature], "model" => ["Trimer"], "list_type" => "LinkedList", "bonds" => create_bond_matrix(N)))
 
 model_matrix = Trimer()
-bonds = create_bond_matrix(N)
-chains = [System(position, species, molecule, density, temperature, model_matrix, bonds; list_type=LinkedList) for _ in 1:M]
-## Define moves and combine them into a pool
+
 pswap = 0.2
 displacement_policy = SimpleGaussian()
 displacement_parameters = ComponentArray(σ=0.05)
 pool = (
-    Move(Displacement(0, zero(box), 0.0), displacement_policy, displacement_parameters, 1.0),
+    Move(Displacement(0, zero(chains[1].box), 0.0), displacement_policy, displacement_parameters, 1.0),
 )
 ## Define the simulation struct
 steps = 1000