Check if initial energy is infinite. Add isinf(e1) in perform_action

Author: romainljsimon

src/atoms.jl

@@ -1,4 +1,4 @@
-struct Atoms{D, Q, V<:AbstractVector, C<:NeighbourList, H, T<:AbstractFloat, SM<:AbstractArray} <: Particles
+struct Atoms{D, V<:AbstractVector, C<:NeighbourList, H, T<:AbstractFloat, SM<:AbstractArray} <: Particles
     position::Vector{SVector{D,T}}
     species::V
     density::T
@@ -8,7 +8,6 @@ struct Atoms{D, Q, V<:AbstractVector, C<:NeighbourList, H, T<:AbstractFloat, SM<
     N::Int
     d::Int
     box::SVector{D,T}
-    local_energy::MVector{Q, T}
     neighbour_list::C
     species_list::H
 end
@@ -20,15 +19,17 @@ function System(position, species, density::T, temperature::T, model_matrix; lis
     d = length(Array(position)[1])
     energy = Vector{T}(undef, 1)
     box = @SVector fill(T((N / density)^(1 / d)), d)
-    local_energy = MVector{N, T}(zeros(T, N))
-    indices = MVector{N, Bool}(falses(N))
     maxcut = maximum([model.rcut for model in model_matrix])
     neighbour_list = list_type(box, maxcut, N)
     species_list = isa(species[1], Integer) ? SpeciesList(species) : nothing
     system = Atoms(position, species, density, energy, temperature, model_matrix, N, d, box, neighbour_list, species_list)
     build_neighbour_list!(system)
-    system.local_energy .= [compute_energy_particle(system, i) for i in eachindex(system)]
-    system.energy[1] = sum(system.local_energy) / 2
+    local_energy = [compute_energy_particle(system, i, neighbour_list) for i in eachindex(position)]
+    energy = mean(local_energy) / 2
+    if isinf(energy) || isnan(energy)
+        error("Initial configuration has infinite or NaN energy.")
+    end
+    system.energy[1] = energy
     return system
 end
 

src/molecules.jl

@@ -13,7 +13,6 @@ struct Molecules{D,  VS<:AbstractVector, C<:NeighbourList, T<:AbstractFloat, SM<
     N::Int
     Nmol::Int
     box::SVector{D,T}
-    local_energy::Vector{T}
     neighbour_list::C
     bonds::Vector{Vector{Int}}
 end
@@ -27,14 +26,17 @@ function System(position, species, molecule, density::T, temperature::T, model_m
     molecule_species = something(molecule_species, ones(Int, N))
     d = length(Array(position)[1])
     box = @SVector fill(T((N / density)^(1 / d)), d)
-    local_energy = zeros(T, N)
     energy = zeros(T, 1)
     maxcut = maximum([model.rcut for model in model_matrix])
     neighbour_list = list_type(box, maxcut, N)
     system = Molecules(position, species, molecule, molecule_species,  start_mol, length_mol, density, temperature, energy, model_matrix, d, N, Nmol,box, neighbour_list, bonds)
     build_neighbour_list!(system)
-    system.local_energy .= [compute_energy_particle(system, i, neighbour_list) for i in eachindex(position)]
-    system.energy[1] = sum(system.local_energy) / 2
+    local_energy = [compute_energy_particle(system, i, neighbour_list) for i in eachindex(position)]
+    energy = mean(local_energy) / 2
+    if isinf(energy) || isnan(energy)
+        error("Initial configuration has infinite or NaN energy.")
+    end
+    system.energy[1] = energy
     return system
 end
 

src/moves.jl

@@ -2,7 +2,7 @@ using ComponentArrays
 
 function Arianna.perform_action!(system::Particles, action::Action)
     e₁, e₂ = perform_action!(system, action)
-    if isinf(e₂)
+    if isinf(e₁) || isinf(e₂)
         action.δe = zero(typeof(system.energy[1]))
     else
         action.δe = e₂ - e₁