Merge pull request #41 from TheDisorderedOrganization/cleanup2

Continuation of the clean up

Author: romainljsimon

src/IO/xyz.jl

@@ -75,7 +75,7 @@ end
 function write_header(io, system::Particles, t, format::XYZ, digits::Integer)
     println(io, length(system))
     box = replace(replace(string(system.box), r"[\[\]]" => ""), r",\s+" => ",")
-    println(io, "step:$t columns:$(get_system_column(system, format))species,position dt:1 cell:$(box) rho:$(system.density) T:$(system.temperature) potential_energy_per_particle:$(mean(system.local_energy)/2)")
+    println(io, "step:$t columns:$(get_system_column(system, format))species,position dt:1 cell:$(box) rho:$(system.density) T:$(system.temperature)")
     return nothing
 end
 

src/ParticlesMC.jl

@@ -16,12 +16,11 @@ include("moves.jl")
 get_position(system::Particles, i::Int) = @inbounds system.position[i]
 get_species(system::Particles, i::Int) = @inbounds system.species[i]
 get_model(system::Particles, i::Int, j::Int) = @inbounds system.model_matrix[get_species(system, i), get_species(system, j)]
-get_local_energy(system::Particles, i::Int) = @inbounds system.local_energy[i]
 get_box(system::Particles) = system.box
 get_neighbour_list(system::Particles) = system.neighbour_list
 Base.length(system::Particles) = system.N
 Base.eachindex(system::Particles) = Base.OneTo(length(system))
-Base.getindex(system::Atoms, i::Int) = system.position[i], system.species[i], system.local_energy[i]
+Base.getindex(system::Atoms, i::Int) = system.position[i], system.species[i]
 
 function Base.iterate(system::Union{Atoms, Molecules}, state=1)
     state > length(system) && return nothing  # Stop iteration

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, local_energy, neighbour_list, species_list)
+    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, local_energy, neighbour_list, bonds)
+    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₁
@@ -141,15 +141,21 @@ end
 struct EnergyBias <: Policy end
 
 function Arianna.log_proposal_density(action::DiscreteSwap, ::EnergyBias, parameters, system::Particles)
-    numerator = parameters.θ₁ * system.local_energy[action.i] + parameters.θ₂ * system.local_energy[action.j]
-    log_sum_exp_1 = log(sum(exp.(parameters.θ₁ .* system.local_energy[system.species_list.sp_ids[action.species[1]]])))
-    log_sum_exp_2 = log(sum(exp.(parameters.θ₂ .* system.local_energy[system.species_list.sp_ids[action.species[2]]])))
+    energy_particle_i = compute_energy_particle(system, action.i, system.neighbour_list)
+    energy_particle_j = compute_energy_particle(system, action.j, system.neighbour_list)
+    numerator = parameters.θ₁ * energy_particle_i+ parameters.θ₂ * energy_particle_j
+    energy_particle_1 = compute_energy_particle(system, system.species_list.sp_ids[action.species[1]], system.neighbour_list)
+    energy_particle_2 = compute_energy_particle(system, system.species_list.sp_ids[action.species[1]], system.neighbour_list)
+    log_sum_exp_1 = log(sum(exp.(parameters.θ₁ .* energy_particle_1)))
+    log_sum_exp_2 = log(sum(exp.(parameters.θ₂ .* energy_particle_1)))
     return numerator - log_sum_exp_1 - log_sum_exp_2
 end
 
 function Arianna.sample_action!(action::DiscreteSwap, ::EnergyBias, parameters, system::Particles, rng)
-    w1s = exp.(parameters.θ₁ .* system.local_energy[system.species_list.sp_ids[action.species[1]]])
-    w2s = exp.(parameters.θ₂ .* system.local_energy[system.species_list.sp_ids[action.species[2]]])
+    energy_particle_1 = compute_energy_particle(system, system.species_list.sp_ids[action.species[1]], system.neighbour_list)
+    energy_particle_2 = compute_energy_particle(system, system.species_list.sp_ids[action.species[1]], system.neighbour_list)
+    w1s = exp.(parameters.θ₁ .* energy_particle_1)
+    w2s = exp.(parameters.θ₂ .* energy_particle_2)
     w1s .= w1s ./ sum(w1s)
     w2s .= w2s ./ sum(w2s)
     id1 = rand(rng, Categorical(w1s))

test/gerhard_energy_distribution.jl

@@ -11,8 +11,8 @@ chains_ll = load_chains("test/config_0.lmp", args=Dict("temperature" => [0.231],
 system = chains[1]
 system_ll = chains_ll[1]
 # GERHARD: -2.676832
-@show mean(system.local_energy) / 2
-@show mean(system_ll.local_energy) / 2
+@show system.energy[1]
+@show system_ll.energy[1]
 chains_bkp = deepcopy(chains)
 chains_ll_bkp = deepcopy(chains_ll)
 

test/runtests.jl

@@ -17,8 +17,8 @@ using DelimitedFiles
     @test system_el.temperature == system_ll.temperature
     @test all.(system_el.position == system_ll.position)
     @test all.(system_el.species == system_ll.species)
-    energy_el = mean(system_el.local_energy) / 2
-    energy_ll = mean(system_ll.local_energy) / 2
+    energy_el = system_el.energy[1]
+    energy_ll = system_ll.energy[1]
     @test isapprox(energy_el, -2.676832, atol=1e-6)
     @test isapprox(energy_ll, -2.676832, atol=1e-6)
 
@@ -118,8 +118,8 @@ end
     @test all.(system_el.position == system_ll.position)
     @test all.(system_el.species == system_ll.species)
     @test all.(system_el.bonds == system_ll.bonds)
-    energy_el = mean(system_el.local_energy) / 2
-    energy_ll = mean(system_ll.local_energy) / 2
+    energy_el = system_el.energy[1]
+    energy_ll = system_ll.energy[1]
     @test isapprox(energy_el, 25.65865662277199, atol=1e-6)
     @test isapprox(energy_ll, 25.65865662277199, atol=1e-6)