Deletion of local_energy variable

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

@@ -25,7 +25,7 @@ function System(position, species, density::T, temperature::T, model_matrix; lis
     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

src/molecules.jl

@@ -31,7 +31,7 @@ function System(position, species, molecule, density::T, temperature::T, model_m
     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

src/moves.jl

@@ -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)