From 6aea28615ae569d247748acf57270b341b7dc095 Mon Sep 17 00:00:00 2001
From: LV-IQ <l236550@dac.unicamp.br>
Date: Fri, 21 Jun 2024 17:57:13 -0300
Subject: [PATCH 01/10] Added min_dist option

---
 src/Reweighting/Reweight.jl    | 38 +++++++++++++++++++++++++---------
 src/Reweighting/Reweighting.jl | 22 +++++++++++++++++++-
 2 files changed, 49 insertions(+), 11 deletions(-)

diff --git a/src/Reweighting/Reweight.jl b/src/Reweighting/Reweight.jl
index 7fe8c7de..cbef841b 100644
--- a/src/Reweighting/Reweight.jl
+++ b/src/Reweighting/Reweight.jl
@@ -138,7 +138,9 @@ function reweight(
     simulation::Simulation, 
     f_perturbation::Function, 
     group_1::AbstractVector{<:Integer}, 
-    group_2::AbstractVector{<:Integer};     
+    group_2::AbstractVector{<:Integer};
+    min_dist::Bool = false,
+    n_atoms_per_molecule::AbstractVector{<:Integer} = [0,0],     
     cutoff::Real = 12.0, 
     k::Real = 1.0, 
     T::Real = 1.0
@@ -150,15 +152,31 @@ function reweight(
         coordinates = positions(frame)
         first_coors = coordinates[group_1]
         second_coors = coordinates[group_2]
-        system = ParticleSystem(
-            xpositions = first_coors,
-            ypositions = second_coors,
-            unitcell = unitcell(frame),
-            cutoff = cutoff,
-            output = 0.0,
-            output_name = :total_energy
-        )
-        energy_vec[iframe] = map_pairwise!((x, y, i, j, d2, total_energy) -> total_energy + f_perturbation(i, j, sqrt(d2)/10), system)
+        system = nothing
+        if min_dist
+            system = minimum_distances(
+                xpositions = coordinates[group_1], 
+                ypositions = coordinates[group_2],
+                unitcell=unitcell(frame), 
+                cutoff = 0.1, 
+                xn_atoms_per_molecule=n_atoms_per_molecule[1],
+                yn_atoms_per_molecule=n_atoms_per_molecule[2]
+            )
+            for mindist in system
+                dist = mindist.d 
+                energy_vec[iframe] += dist -> f_perturbation(dist)
+            end
+        else
+            system = ParticleSystem(
+                xpositions = first_coors,
+                ypositions = second_coors,
+                unitcell = unitcell(frame),
+                cutoff = cutoff,
+                output = 0.0,
+                output_name = :total_energy
+            )
+            energy_vec[iframe] = map_pairwise!((x, y, i, j, d2, total_energy) -> total_energy + f_perturbation(i, j, sqrt(d2)/10), system)
+        end
     end
     @. prob_rel_vec = exp(-(energy_vec)/k*T)
     prob_vec = prob_rel_vec/sum(prob_rel_vec)
diff --git a/src/Reweighting/Reweighting.jl b/src/Reweighting/Reweighting.jl
index 7936f252..ac99714d 100644
--- a/src/Reweighting/Reweighting.jl
+++ b/src/Reweighting/Reweighting.jl
@@ -1,7 +1,7 @@
 module Reweighting
 
 using CellListMap: ParticleSystem, map_pairwise, map_pairwise!
-using ..MolSimToolkit: Simulation, positions, unitcell
+using ..MolSimToolkit: Simulation, positions, unitcell, minimum_distances
 
 export reweight, lennard_jones_perturbation, poly_decay_perturbation, gaussian_decay_perturbation
 
@@ -77,4 +77,24 @@ end
         0.09973413264337352
         0.12988266765019355
     ]
+end
+
+@testitem "Reweighting small trajectory using minimum distances" begin
+    import PDBTools
+    using MolSimToolkit.Reweighting
+    using MolSimToolkit.Reweighting: testdir
+
+    simulation = Simulation("$testdir/Testing_reweighting.pdb", "$testdir/Testing_reweighting_10_frames_trajectory.xtc")
+
+    i1 = PDBTools.selindex(atoms(simulation), "index 97 or index 106")
+
+    i2 = PDBTools.selindex(atoms(simulation), "residue 15 and name HB3")
+
+    sum_of_dist = reweight(simulation, (r) -> r, i1, i2; min_dist = true, n_atoms_per_molecule = [2,10], cutoff = 25.0)
+    @test sum_of_dist.energy ≈ [
+        1.773896547670759, 1.5923698293115915, 1.716614676290554, 
+        1.933003841107648, 1.602329229247863, 1.9639005665480983, 
+        3.573986006775934, 2.188798265022823, 2.066180657974777, 
+        1.6845109623700647
+    ]
 end
\ No newline at end of file

From 5a72a7c8452a6524f3560583dadb0eaf1dbdcd96 Mon Sep 17 00:00:00 2001
From: LV-IQ <l236550@dac.unicamp.br>
Date: Mon, 1 Jul 2024 09:14:03 -0300
Subject: [PATCH 02/10] Fixed file functions

---
 src/Reweighting/Reweight.jl    | 22 +++++++++++-----------
 src/Reweighting/Reweighting.jl | 23 +++++++++++++++--------
 2 files changed, 26 insertions(+), 19 deletions(-)

diff --git a/src/Reweighting/Reweight.jl b/src/Reweighting/Reweight.jl
index cbef841b..31c7fab6 100644
--- a/src/Reweighting/Reweight.jl
+++ b/src/Reweighting/Reweight.jl
@@ -140,7 +140,7 @@ function reweight(
     group_1::AbstractVector{<:Integer}, 
     group_2::AbstractVector{<:Integer};
     min_dist::Bool = false,
-    n_atoms_per_molecule::AbstractVector{<:Integer} = [0,0],     
+    n_atoms_per_molecule::Real = 0,     
     cutoff::Real = 12.0, 
     k::Real = 1.0, 
     T::Real = 1.0
@@ -152,19 +152,19 @@ function reweight(
         coordinates = positions(frame)
         first_coors = coordinates[group_1]
         second_coors = coordinates[group_2]
-        system = nothing
         if min_dist
-            system = minimum_distances(
-                xpositions = coordinates[group_1], 
-                ypositions = coordinates[group_2],
+            minimum_dists = minimum_distances(
+                xpositions = first_coors, 
+                ypositions = second_coors,
                 unitcell=unitcell(frame), 
-                cutoff = 0.1, 
-                xn_atoms_per_molecule=n_atoms_per_molecule[1],
-                yn_atoms_per_molecule=n_atoms_per_molecule[2]
+                cutoff = cutoff, 
+                xn_atoms_per_molecule=n_atoms_per_molecule,
             )
-            for mindist in system
-                dist = mindist.d 
-                energy_vec[iframe] += dist -> f_perturbation(dist)
+            for mindist in minimum_dists
+                if mindist.within_cutoff
+                    dist = mindist.d 
+                    energy_vec[iframe] += f_perturbation(dist)
+                end
             end
         else
             system = ParticleSystem(
diff --git a/src/Reweighting/Reweighting.jl b/src/Reweighting/Reweighting.jl
index ac99714d..8aee09b6 100644
--- a/src/Reweighting/Reweighting.jl
+++ b/src/Reweighting/Reweighting.jl
@@ -1,7 +1,8 @@
 module Reweighting
 
 using CellListMap: ParticleSystem, map_pairwise, map_pairwise!
-using ..MolSimToolkit: Simulation, positions, unitcell, minimum_distances
+using ..MolSimToolkit: Simulation, positions, unitcell
+using ..MolSimToolkit.MolecularMinimumDistances
 
 export reweight, lennard_jones_perturbation, poly_decay_perturbation, gaussian_decay_perturbation
 
@@ -86,15 +87,21 @@ end
 
     simulation = Simulation("$testdir/Testing_reweighting.pdb", "$testdir/Testing_reweighting_10_frames_trajectory.xtc")
 
-    i1 = PDBTools.selindex(atoms(simulation), "index 97 or index 106")
+    i1 = PDBTools.selindex(atoms(simulation), "resname TFE and name O")
 
-    i2 = PDBTools.selindex(atoms(simulation), "residue 15 and name HB3")
+    i2 = PDBTools.selindex(atoms(simulation), "residue 7")
 
-    sum_of_dist = reweight(simulation, (r) -> r, i1, i2; min_dist = true, n_atoms_per_molecule = [2,10], cutoff = 25.0)
+    sum_of_dist = reweight(simulation, (r) -> r, i1, i2; min_dist = true, n_atoms_per_molecule = 7, cutoff = 10.0)
     @test sum_of_dist.energy ≈ [
-        1.773896547670759, 1.5923698293115915, 1.716614676290554, 
-        1.933003841107648, 1.602329229247863, 1.9639005665480983, 
-        3.573986006775934, 2.188798265022823, 2.066180657974777, 
-        1.6845109623700647
+    62.492586115838314, 
+    74.1608164984956, 
+    74.7700517897278, 
+    50.851563012130356, 
+    62.994684516904, 
+    67.27903007182715, 
+    81.47993128439455, 
+    52.87729020037196, 
+    63.76631345255868, 
+    55.55921195504992
     ]
 end
\ No newline at end of file

From df65be4ae5f819dd297340f14893ba4a3f5866ad Mon Sep 17 00:00:00 2001
From: LV-IQ <l236550@dac.unicamp.br>
Date: Mon, 1 Jul 2024 09:57:26 -0300
Subject: [PATCH 03/10] Added mindist for one set

---
 src/Reweighting/Reweight.jl    | 41 +++++++++++++++++++++---------
 src/Reweighting/Reweighting.jl | 46 ++++++++++++++++++++++++++--------
 2 files changed, 64 insertions(+), 23 deletions(-)

diff --git a/src/Reweighting/Reweight.jl b/src/Reweighting/Reweight.jl
index 31c7fab6..372417c5 100644
--- a/src/Reweighting/Reweight.jl
+++ b/src/Reweighting/Reweight.jl
@@ -109,7 +109,9 @@ This result is the energy difference between the  perturbed frame and the origin
 function reweight(
     simulation::Simulation, 
     f_perturbation::Function, 
-    group_1::AbstractVector{<:Integer}; 
+    group_1::AbstractVector{<:Integer};
+    min_dist::Bool = false,
+    n_atoms_per_molecule::Real = 0, 
     cutoff::Real = 12.0, 
     k::Real = 1.0, 
     T::Real = 1.0
@@ -120,14 +122,29 @@ function reweight(
     for (iframe, frame) in enumerate(simulation)
         coordinates = positions(frame)
         first_coors = coordinates[group_1]
-        system = ParticleSystem(
-            xpositions = first_coors,
-            unitcell = unitcell(frame),
-            cutoff = cutoff,
-            output = 0.0,
-            output_name = :total_energy
-        )
-        energy_vec[iframe] = map_pairwise!((x, y, i, j, d2, total_energy) -> total_energy + f_perturbation(i, j, sqrt(d2)/10), system)
+        if min_dist
+            minimum_dists = SelfPairs(
+                xpositions = first_coors,
+                cutoff = cutoff,
+                unitcell = unitcell(frame),
+                xn_atoms_per_molecule = n_atoms_per_molecule
+            )
+            for mindist in minimum_distances!(minimum_dists)
+                if mindist.within_cutoff
+                    dist = mindist.d 
+                    energy_vec[iframe] += f_perturbation(dist)
+                end
+            end
+        else
+            system = ParticleSystem(
+                xpositions = first_coors,
+                unitcell = unitcell(frame),
+                cutoff = cutoff,
+                output = 0.0,
+                output_name = :total_energy
+            )
+            energy_vec[iframe] = map_pairwise!((x, y, i, j, d2, total_energy) -> total_energy + f_perturbation(i, j, sqrt(d2)/10), system)
+        end
     end
     @. prob_rel_vec = exp(-(energy_vec)/k*T)
     prob_vec = prob_rel_vec/sum(prob_rel_vec)
@@ -156,10 +173,10 @@ function reweight(
             minimum_dists = minimum_distances(
                 xpositions = first_coors, 
                 ypositions = second_coors,
-                unitcell=unitcell(frame), 
+                unitcell = unitcell(frame), 
                 cutoff = cutoff, 
-                xn_atoms_per_molecule=n_atoms_per_molecule,
-            )
+                xn_atoms_per_molecule = n_atoms_per_molecule,
+            ) 
             for mindist in minimum_dists
                 if mindist.within_cutoff
                     dist = mindist.d 
diff --git a/src/Reweighting/Reweighting.jl b/src/Reweighting/Reweighting.jl
index 8aee09b6..0168a58c 100644
--- a/src/Reweighting/Reweighting.jl
+++ b/src/Reweighting/Reweighting.jl
@@ -91,17 +91,41 @@ end
 
     i2 = PDBTools.selindex(atoms(simulation), "residue 7")
 
-    sum_of_dist = reweight(simulation, (r) -> r, i1, i2; min_dist = true, n_atoms_per_molecule = 7, cutoff = 10.0)
+    sum_of_dist = reweight(simulation, (r) -> r, i1, i2; min_dist = true, n_atoms_per_molecule = 1, cutoff = 10.0)
     @test sum_of_dist.energy ≈ [
-    62.492586115838314, 
-    74.1608164984956, 
-    74.7700517897278, 
-    50.851563012130356, 
-    62.994684516904, 
-    67.27903007182715, 
-    81.47993128439455, 
-    52.87729020037196, 
-    63.76631345255868, 
-    55.55921195504992
+        81.63048517945049, 
+        81.75153397197373, 
+        83.4725848687559, 
+        50.851563012130356, 
+        70.51272314980994, 
+        67.27903007182715, 
+        81.47993128439455, 
+        71.87915113998032, 
+        72.327486822439, 
+        55.55921195504992
+    ]
+end
+
+@testitem "Reweighting small trajectory using minimum distances with one set" begin
+    import PDBTools
+    using MolSimToolkit.Reweighting
+    using MolSimToolkit.Reweighting: testdir
+
+    simulation = Simulation("$testdir/Testing_reweighting.pdb", "$testdir/Testing_reweighting_10_frames_trajectory.xtc")
+
+    i1 = PDBTools.selindex(atoms(simulation), "resname TFE and name O")
+
+    sum_of_dist = reweight(simulation, (r) -> r, i1; min_dist = true, n_atoms_per_molecule = 1, cutoff = 10.0)
+    @test sum_of_dist.energy ≈ [
+        1461.0357719672927, 
+        1441.135120247483, 
+        1459.9465255299242, 
+        1475.6216777213508, 
+        1474.9315243351775, 
+        1474.7965527585634, 
+        1452.6348409088291, 
+        1478.65551753013, 
+        1452.50720582719, 
+        1507.8298874600328
     ]
 end
\ No newline at end of file

From 07b5d11af042f264dd4026372f0e51416ed92b78 Mon Sep 17 00:00:00 2001
From: CompatHelper Julia <compathelper_noreply@julialang.org>
Date: Mon, 8 Jul 2024 01:16:09 +0000
Subject: [PATCH 04/10] CompatHelper: bump compat for TestItems to 1, (keep
 existing compat)

---
 Project.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Project.toml b/Project.toml
index b4aa9a61..3d827cfe 100644
--- a/Project.toml
+++ b/Project.toml
@@ -51,7 +51,7 @@ Statistics = "1.9"
 StatsBase = "0.34"
 Test = "1.9"
 TestItemRunner = "0.2.1"
-TestItems = "0.1.1"
+TestItems = "0.1.1, 1"
 julia = "1.9"
 
 [extras]

From f9ce419da31546ed2f96e709b4a6c3932e50bb16 Mon Sep 17 00:00:00 2001
From: LV-IQ <l236550@dac.unicamp.br>
Date: Tue, 23 Jul 2024 16:06:51 -0300
Subject: [PATCH 05/10] Text and version fix (not finished)

---
 Project.toml            |  4 ++--
 docs/make.jl            |  2 +-
 docs/src/Reweighting.md | 34 +++++++++++++++++++++++-----------
 3 files changed, 26 insertions(+), 14 deletions(-)

diff --git a/Project.toml b/Project.toml
index 3d827cfe..f0e7762a 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "MolSimToolkit"
 uuid = "054db54f-6694-444d-9bbb-e9ecdbfe77be"
 authors = ["Leandro Martinez <lmartine@unicamp.br> and contributors"]
-version = "1.14.1-DEV"
+version = "1.16.1-DEV"
 
 [deps]
 AtomsBase = "a963bdd2-2df7-4f54-a1ee-49d51e6be12a"
@@ -31,7 +31,7 @@ Plotting = "Plots"
 Aqua = "0.8"
 AtomsBase = "0.3.4"
 BenchmarkTools = "1.3"
-CellListMap = "0.8.30, 0.9"
+CellListMap = "0.9"
 Chemfiles = "0.10.31"
 DelimitedFiles = "1.9"
 DocStringExtensions = "0.9"
diff --git a/docs/make.jl b/docs/make.jl
index d7b08b07..6ad1a9d9 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -19,7 +19,7 @@ makedocs(
         "Miscelaneous functions" => "miscelaneous.md",
         "Developer zone" => "Developer.md",
         "Plotting style" => "plotting_style.md",
-        "Simulation Reweight" => "Reweighting.md"
+        "Simulation Reweighting" => "Reweighting.md"
     ],
 )
 deploydocs(
diff --git a/docs/src/Reweighting.md b/docs/src/Reweighting.md
index 3da17dd4..3cc64b81 100644
--- a/docs/src/Reweighting.md
+++ b/docs/src/Reweighting.md
@@ -34,40 +34,47 @@ julia> i1 = PDBTools.selindex(atoms(simulation), "resname TFE and name O")
 
 julia> i2 = PDBTools.selindex(atoms(simulation), "protein and name O")
 ```
+In the picture below, it is shown the interaction which will be perturbed in this system:
+
+!!!!INSERT FIGURE
 
 ## Setting perturbation function
 In order to obtain these weights, we have to use two functions: the ```reweight``` function, which will calculate each weight and the ```perturbation``` function, responsible for taking each computed distance between atomic pairs in every frame and determine the resulting energy using these distances in that particular frame based on the applied perturbation.
 
-So, secondly, we define some "perturbation" function (here we call it ```gaussian decay```) and set up its parameters. Please, take a look at the interface:
+So, secondly, we define some "perturbation" function (here we call it ```polynomial_decay```) and set up its parameters. Please, take a look at the interface:
 
 ```julia-repl
-julia> gaussian_decay(r, α, β) = α*exp(-abs(β)*r^2)
+julia> poly_decay_perturbation(r, α, cut) = α * ((r/cut)^2 - 1)^2
 gaussian_decay (generic function with 1 method)
 
 julia> α = 5.e-3
 0.005
 
-julia> β = 5.e-3
-0.005
+julia> cut = 10.0
+10.0
 ```
 
-As it can be seen, the function has to receive two parameters: `r` which corresponds to the distance between two selected atoms and some parameter to account a modification and change its magnitude, here, we inserted two of them in the same function, `α`, to change the maximum value of the gaussian curve and `β`, to adjust its decay behaviour with a given value of `r`.
+
+As it can be seen, the function has to receive three parameters: `r` which corresponds to the distance between two selected atoms and some parameter to account a modification and change its magnitude, here, we inserted two of them in the same function, `α`, to change the maximum value of the curve (at r = 1) and `cut`, the distance `r` where the function equals zero. In the image below, we can see the curve and how it changes with different values of `α` and `cut`:
+
+!!! INSERT FIGURE
 
 ## Computing the new weights
 Finally, using the ```reweight``` function, we pass both the ```simulation``` and the last function anonymously in the input. Again, watch the interface:
 
 ```julia-repl
-julia> cut_off = 12.0
-12.0
 
-julia> weights = reweight(simulation, (i,j,r) -> gaussian_decay(r, α, β), i1, i2; cutoff = cut_off)
+julia> weights = reweight(simulation, r -> gaussian_decay(r, α, cut), i1, i2; cutoff = cut)
 ```
 
-`i and j`: if you selected two atom types, `i` will be the index for either the first, the second, the third and so on up to the last atom of the first group and `j` will be same, but now for the second one. With these two parameters, it is possible to determine every combination of two atoms, each one coming from one group, and compute the associated dsitance `r`, so that we are taking into account all interactions between these two atom types to our perturbation. However, if we are dealing with just one group, both of them are indexes for all the atoms of the selected group. Bear in mind that repeated combinations (like `i,j = 1,2 or 2,1`) will no be computed, since the `reweight` function calls the `map_pairwise!` function, from [CellListMap.jl](https://github.com/m3g/CellListMap.jl) that is able to avoid this problem.
+`r`: the distance between the twos atoms 
 
-`r`: the distance between the twos atoms with indexes `i` and `j` in the selected groups.
+`cutoff`: the maximum distance that will be computed between two atoms. The default value is `12.0 Å`.
 
-`cutoff`: the maximum distance that will be computed between two atoms. The default value is `12.0` Angstrom.
+!!! warning
+    It is highly recommended to set the same value of `cutoff` for both `perturbation` and `reweight` functions.
+    With this in mind, calculations will be done more quickly and you do not need to worry about your input function
+    behaviour above the `cutoff` value, since distances out of the perturbation range will not be computed.
 
 Once the calculations are finished, the resulted interface is shown, like the example below:
 ```julia-repl
@@ -120,6 +127,11 @@ julia> weights.probability
  0.12988266765019355
 ```
 
+!!! tip
+    Note that these values (and, consequently, calculations that use them) are functions of `r`, 
+    so in other to avoid mathematical complications, a good piece of advice is to create functions that 
+    are continous in the closed interval from zero to the cutoff value.
+
 ## Reference Functions
 ```@autodocs
 Modules = [MolSimToolkit.Reweighting]

From 6fe7976640690cb2e74401313b961991e8bf4e68 Mon Sep 17 00:00:00 2001
From: LV-IQ <l236550@dac.unicamp.br>
Date: Tue, 23 Jul 2024 16:07:20 -0300
Subject: [PATCH 06/10] Test examples fix

---
 src/Reweighting/Reweight.jl              | 79 +++++++-----------------
 src/Reweighting/Reweighting.jl           | 70 ++++-----------------
 src/Reweighting/perturbation_examples.jl | 26 ++++----
 3 files changed, 47 insertions(+), 128 deletions(-)

diff --git a/src/Reweighting/Reweight.jl b/src/Reweighting/Reweight.jl
index 372417c5..988fdf5c 100644
--- a/src/Reweighting/Reweight.jl
+++ b/src/Reweighting/Reweight.jl
@@ -68,7 +68,7 @@ julia> i2 = PDBTools.selindex(atoms(simulation), "residue 15 and name HB3")
 1-element AbstractVector{<:Integer}:
  171
 
-julia> sum_of_dist = reweight(simulation, (i,j,d2) -> d2, i1, i2; cutoff = 25.0)
+julia> sum_of_dist = reweight(simulation, r -> r, i1, i2; cutoff = 25.0)
 -------------
 FRAME WEIGHTS
 -------------
@@ -110,8 +110,6 @@ function reweight(
     simulation::Simulation, 
     f_perturbation::Function, 
     group_1::AbstractVector{<:Integer};
-    min_dist::Bool = false,
-    n_atoms_per_molecule::Real = 0, 
     cutoff::Real = 12.0, 
     k::Real = 1.0, 
     T::Real = 1.0
@@ -122,31 +120,16 @@ function reweight(
     for (iframe, frame) in enumerate(simulation)
         coordinates = positions(frame)
         first_coors = coordinates[group_1]
-        if min_dist
-            minimum_dists = SelfPairs(
-                xpositions = first_coors,
-                cutoff = cutoff,
-                unitcell = unitcell(frame),
-                xn_atoms_per_molecule = n_atoms_per_molecule
-            )
-            for mindist in minimum_distances!(minimum_dists)
-                if mindist.within_cutoff
-                    dist = mindist.d 
-                    energy_vec[iframe] += f_perturbation(dist)
-                end
-            end
-        else
-            system = ParticleSystem(
-                xpositions = first_coors,
-                unitcell = unitcell(frame),
-                cutoff = cutoff,
-                output = 0.0,
-                output_name = :total_energy
-            )
-            energy_vec[iframe] = map_pairwise!((x, y, i, j, d2, total_energy) -> total_energy + f_perturbation(i, j, sqrt(d2)/10), system)
-        end
+        system = ParticleSystem(
+            xpositions = first_coors,
+            unitcell = unitcell(frame),
+            cutoff = cutoff,
+            output = 0.0,
+            output_name = :total_energy
+        )
+        energy_vec[iframe] = map_pairwise!((x, y, i, j, d2, total_energy) -> total_energy + f_perturbation(sqrt(d2)), system)
     end
-    @. prob_rel_vec = exp(-(energy_vec)/k*T)
+    @. prob_rel_vec = exp(-energy_vec/(k*T))
     prob_vec = prob_rel_vec/sum(prob_rel_vec)
     output = ReweightResults(prob_vec, prob_rel_vec, energy_vec)
     return output
@@ -155,11 +138,9 @@ function reweight(
     simulation::Simulation, 
     f_perturbation::Function, 
     group_1::AbstractVector{<:Integer}, 
-    group_2::AbstractVector{<:Integer};
-    min_dist::Bool = false,
-    n_atoms_per_molecule::Real = 0,     
+    group_2::AbstractVector{<:Integer};    
     cutoff::Real = 12.0, 
-    k::Real = 1.0, 
+    k::Real = 1.0,
     T::Real = 1.0
 )
     prob_vec = zeros(length(simulation))
@@ -169,33 +150,17 @@ function reweight(
         coordinates = positions(frame)
         first_coors = coordinates[group_1]
         second_coors = coordinates[group_2]
-        if min_dist
-            minimum_dists = minimum_distances(
-                xpositions = first_coors, 
-                ypositions = second_coors,
-                unitcell = unitcell(frame), 
-                cutoff = cutoff, 
-                xn_atoms_per_molecule = n_atoms_per_molecule,
-            ) 
-            for mindist in minimum_dists
-                if mindist.within_cutoff
-                    dist = mindist.d 
-                    energy_vec[iframe] += f_perturbation(dist)
-                end
-            end
-        else
-            system = ParticleSystem(
-                xpositions = first_coors,
-                ypositions = second_coors,
-                unitcell = unitcell(frame),
-                cutoff = cutoff,
-                output = 0.0,
-                output_name = :total_energy
-            )
-            energy_vec[iframe] = map_pairwise!((x, y, i, j, d2, total_energy) -> total_energy + f_perturbation(i, j, sqrt(d2)/10), system)
-        end
+        system = ParticleSystem(
+            xpositions = first_coors,
+            ypositions = second_coors,
+            unitcell = unitcell(frame),
+            cutoff = cutoff,
+            output = 0.0,
+            output_name = :total_energy
+        )
+        energy_vec[iframe] = map_pairwise!((x, y, i, j, d2, total_energy) -> total_energy + f_perturbation(sqrt(d2)), system)
     end
-    @. prob_rel_vec = exp(-(energy_vec)/k*T)
+    @. prob_rel_vec = exp(-energy_vec/(k*T))
     prob_vec = prob_rel_vec/sum(prob_rel_vec)
     output = ReweightResults(prob_vec, prob_rel_vec, energy_vec)
     return output
diff --git a/src/Reweighting/Reweighting.jl b/src/Reweighting/Reweighting.jl
index 0168a58c..5bddfcb1 100644
--- a/src/Reweighting/Reweighting.jl
+++ b/src/Reweighting/Reweighting.jl
@@ -24,7 +24,7 @@ end #Module Reweighting
 
     i2 = PDBTools.selindex(atoms(simulation), "residue 11")
 
-    sum_of_dist = reweight(simulation, (i,j,r) -> r, [i1[239]], i2; cutoff = 25.0)
+    sum_of_dist = reweight(simulation, r -> r/10, [i1[239]], i2; cutoff = 25.0)
     @test sum_of_dist.energy ≈ [7.4295543149]
 end
 
@@ -39,11 +39,17 @@ end
 
     i2 = PDBTools.selindex(atoms(simulation), "residue 15 and name HB3")
 
-    sum_of_dist = reweight(simulation, (i,j,r) -> r, i1, i2, cutoff = 25.0)
+    sum_of_dist = reweight(simulation, r -> r/10, i1, i2, cutoff = 25.0)
     @test sum_of_dist.energy ≈ [
-        1.773896547670759, 1.5923698293115915, 1.716614676290554, 
-        1.933003841107648, 1.602329229247863, 1.9639005665480983, 
-        3.573986006775934, 2.188798265022823, 2.066180657974777, 
+        1.773896547670759, 
+        1.5923698293115915, 
+        1.716614676290554, 
+        1.933003841107648, 
+        1.602329229247863, 
+        1.9639005665480983, 
+        3.573986006775934, 
+        2.188798265022823, 
+        2.066180657974777, 
         1.6845109623700647
     ]
 end
@@ -63,9 +69,7 @@ end
 
     β = 5.e-3
 
-    cut_off = 12.0
-
-    probs_test = reweight(simulation, (i,j,r) -> gaussian_decay_perturbation(r, α, β), i1, i2; cutoff = cut_off)
+    probs_test = reweight(simulation, r -> gaussian_decay_perturbation(r/10, α, β), i1, i2; cutoff = 12.0)
     @test probs_test.probability ≈ [
         0.08987791339898044
         0.07326337222373071
@@ -78,54 +82,4 @@ end
         0.09973413264337352
         0.12988266765019355
     ]
-end
-
-@testitem "Reweighting small trajectory using minimum distances" begin
-    import PDBTools
-    using MolSimToolkit.Reweighting
-    using MolSimToolkit.Reweighting: testdir
-
-    simulation = Simulation("$testdir/Testing_reweighting.pdb", "$testdir/Testing_reweighting_10_frames_trajectory.xtc")
-
-    i1 = PDBTools.selindex(atoms(simulation), "resname TFE and name O")
-
-    i2 = PDBTools.selindex(atoms(simulation), "residue 7")
-
-    sum_of_dist = reweight(simulation, (r) -> r, i1, i2; min_dist = true, n_atoms_per_molecule = 1, cutoff = 10.0)
-    @test sum_of_dist.energy ≈ [
-        81.63048517945049, 
-        81.75153397197373, 
-        83.4725848687559, 
-        50.851563012130356, 
-        70.51272314980994, 
-        67.27903007182715, 
-        81.47993128439455, 
-        71.87915113998032, 
-        72.327486822439, 
-        55.55921195504992
-    ]
-end
-
-@testitem "Reweighting small trajectory using minimum distances with one set" begin
-    import PDBTools
-    using MolSimToolkit.Reweighting
-    using MolSimToolkit.Reweighting: testdir
-
-    simulation = Simulation("$testdir/Testing_reweighting.pdb", "$testdir/Testing_reweighting_10_frames_trajectory.xtc")
-
-    i1 = PDBTools.selindex(atoms(simulation), "resname TFE and name O")
-
-    sum_of_dist = reweight(simulation, (r) -> r, i1; min_dist = true, n_atoms_per_molecule = 1, cutoff = 10.0)
-    @test sum_of_dist.energy ≈ [
-        1461.0357719672927, 
-        1441.135120247483, 
-        1459.9465255299242, 
-        1475.6216777213508, 
-        1474.9315243351775, 
-        1474.7965527585634, 
-        1452.6348409088291, 
-        1478.65551753013, 
-        1452.50720582719, 
-        1507.8298874600328
-    ]
 end
\ No newline at end of file
diff --git a/src/Reweighting/perturbation_examples.jl b/src/Reweighting/perturbation_examples.jl
index c11ba23d..409826ca 100644
--- a/src/Reweighting/perturbation_examples.jl
+++ b/src/Reweighting/perturbation_examples.jl
@@ -1,4 +1,4 @@
-#SETTING EQUATIONS FOR PERTUBATION
+#FUNCTIONS FOR PERTURBATIONS
 
 #Perturbing Lennard-Jones interactions
 function switching_func(r, switch, cut)
@@ -6,27 +6,27 @@ function switching_func(r, switch, cut)
     return switching
 end
 
-function lennard_jones(r, sig, ep) 
-    V = 4 * ep * ((sig/r)^12 - (sig/r)^6)
+function lennard_jones(r, σ, ϵ) 
+    V = 4 * ϵ * ((σ/r)^12 - (σ/r)^6)
     return V
 end
 
-function lennard_jones_perturbation(r, sig1, ep1, sig2, ep2; alpha = 0.0, beta = 0.0, smooth::Bool = false, switch = 8.0, cut = 12.0)
-    epsilon_original = sqrt(ep1*ep2)
-    epsilon_perturbed = sqrt(ep1*(1+beta)*ep2)
-    sigma_original = (sig1 + sig2)/2
-    sigma_perturbed = (sig1*(1 + alpha) + sig2)/2
-    V_perturbed = (lennard_jones(r, sigma_perturbed,  epsilon_perturbed) - lennard_jones(r,  sigma_original, epsilon_original))
+function lennard_jones_perturbation(r, σ1, ϵ1, σ2, ϵ2; α = 0.0, β = 0.0, smooth::Bool = false, switch = 8.0, cut = 12.0)
+    ϵ_orig = sqrt(ϵ1*ϵ2)
+    ϵ_pert = (1 + β) * sqrt(ϵ1*ϵ2)
+    σ_orig = (σ1 + σ2)/2
+    σ_pert = α * (σ1* + σ2)/2
+    V_pert = (lennard_jones(r, σ_pert,  ϵ_pert) - lennard_jones(r,  σ_orig, ϵ_orig))
     if smooth && switch < r <= cut
-        V_perturbed = switching_func(r,switch,cut) * V_perturbed
+        V_pert = switching_func(r, switch, cut) * V_pert
     else
-        V_perturbed =  V_perturbed - (lennard_jones(cut, sigma_perturbed,  epsilon_perturbed) - lennard_jones(cut,  sigma_original, epsilon_original))
+        V_pert =  V_pert - (lennard_jones(cut, σ_pert,  ϵ_pert) - lennard_jones(cut,  σ_orig, ϵ_orig))
     end
     return V_perturbed
 end
 
 #Applying a polynomial decay pertubation
-poly_decay_perturbation(r, alpha, cut) = r > cut ? zero(r) : alpha*((r/cut)^2 - 1)^2
+poly_decay_perturbation(r, α, cut) = r > cut ? zero(r) : α * ((r/cut)^2 - 1)^2
 
 #Applying a half-gaussian pertubation
-gaussian_decay_perturbation(r, alpha, beta) = beta > 0 ? alpha*exp(-beta*r^2) : error("you need to insert positive values for β")
\ No newline at end of file
+gaussian_decay_perturbation(r, alpha, beta) = beta > 0 ? alpha*exp(-beta*r^2) : error("you need to insert positive values for β")

From 2a50729d792d47711400339482aeb6692db09602 Mon Sep 17 00:00:00 2001
From: LV-IQ <l236550@dac.unicamp.br>
Date: Tue, 23 Jul 2024 16:32:33 -0300
Subject: [PATCH 07/10] Text testing

---
 docs/src/Reweighting.md | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/docs/src/Reweighting.md b/docs/src/Reweighting.md
index c6432b2a..ee57f608 100644
--- a/docs/src/Reweighting.md
+++ b/docs/src/Reweighting.md
@@ -76,7 +76,7 @@ julia> weights = reweight(simulation, r -> gaussian_decay(r, α, cut), i1, i2; c
 
 `cutoff`: the maximum distance that will be computed between two atoms. The default value is `12.0 Å`.
 
-!!! warning
+!!!warning
     It is highly recommended to set the same value of `cutoff` for both `perturbation` and `reweight` functions.
     With this in mind, calculations will be done more quickly and you do not need to worry about your input function
     behaviour above the `cutoff` value, since distances out of the perturbation range will not be computed.
@@ -132,7 +132,7 @@ julia> weights.probability
  0.12988266765019355
 ```
 
-!!! tip
+!!!tip
     Note that these values (and, consequently, calculations that use them) are functions of `r`, 
     so in other to avoid mathematical complications, a good piece of advice is to create functions that 
     are continous in the closed interval from zero to the cutoff value.

From 55c12aa0805545380431ea8ccb4e0926a0cd91de Mon Sep 17 00:00:00 2001
From: LV-IQ <l236550@dac.unicamp.br>
Date: Tue, 23 Jul 2024 16:39:18 -0300
Subject: [PATCH 08/10] Text fixing

---
 docs/src/Reweighting.md | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/docs/src/Reweighting.md b/docs/src/Reweighting.md
index ee57f608..a65f9493 100644
--- a/docs/src/Reweighting.md
+++ b/docs/src/Reweighting.md
@@ -76,7 +76,7 @@ julia> weights = reweight(simulation, r -> gaussian_decay(r, α, cut), i1, i2; c
 
 `cutoff`: the maximum distance that will be computed between two atoms. The default value is `12.0 Å`.
 
-!!!warning
+!!! warning
     It is highly recommended to set the same value of `cutoff` for both `perturbation` and `reweight` functions.
     With this in mind, calculations will be done more quickly and you do not need to worry about your input function
     behaviour above the `cutoff` value, since distances out of the perturbation range will not be computed.
@@ -132,10 +132,10 @@ julia> weights.probability
  0.12988266765019355
 ```
 
-!!!tip
+!!! tip
     Note that these values (and, consequently, calculations that use them) are functions of `r`, 
     so in other to avoid mathematical complications, a good piece of advice is to create functions that 
-    are continous in the closed interval from zero to the cutoff value.
+    are continous in the closed interval from zero to the cutoff value.  
 
 ## Reference Functions
 ```@autodocs

From b4e1df3cd170f524418364b9ce04af16d7de0db5 Mon Sep 17 00:00:00 2001
From: LV-IQ <l236550@dac.unicamp.br>
Date: Fri, 9 Aug 2024 12:25:37 -0300
Subject: [PATCH 09/10] function fix

---
 src/Reweighting/perturbation_examples.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/Reweighting/perturbation_examples.jl b/src/Reweighting/perturbation_examples.jl
index 409826ca..6846ff50 100644
--- a/src/Reweighting/perturbation_examples.jl
+++ b/src/Reweighting/perturbation_examples.jl
@@ -22,7 +22,7 @@ function lennard_jones_perturbation(r, σ1, ϵ1, σ2, ϵ2; α = 0.0, β = 0.0, s
     else
         V_pert =  V_pert - (lennard_jones(cut, σ_pert,  ϵ_pert) - lennard_jones(cut,  σ_orig, ϵ_orig))
     end
-    return V_perturbed
+    return V_pert
 end
 
 #Applying a polynomial decay pertubation

From b439ccd777d8a053f2a26e1f01ee1713d5d3f4ba Mon Sep 17 00:00:00 2001
From: LV-IQ <l236550@dac.unicamp.br>
Date: Tue, 17 Sep 2024 21:55:43 -0300
Subject: [PATCH 10/10] Fixed version

---
 Project.toml            | 2 +-
 docs/src/Reweighting.md | 1 -
 2 files changed, 1 insertion(+), 2 deletions(-)

diff --git a/Project.toml b/Project.toml
index f0e7762a..abd48ff8 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "MolSimToolkit"
 uuid = "054db54f-6694-444d-9bbb-e9ecdbfe77be"
 authors = ["Leandro Martinez <lmartine@unicamp.br> and contributors"]
-version = "1.16.1-DEV"
+version = "1.19.1-DEV"
 
 [deps]
 AtomsBase = "a963bdd2-2df7-4f54-a1ee-49d51e6be12a"
diff --git a/docs/src/Reweighting.md b/docs/src/Reweighting.md
index a65f9493..94354694 100644
--- a/docs/src/Reweighting.md
+++ b/docs/src/Reweighting.md
@@ -59,7 +59,6 @@ julia> cut = 10.0
 10.0
 ```
 
-
 As it can be seen, the function has to receive three parameters: `r` which corresponds to the distance between two selected atoms and some parameter to account a modification and change its magnitude, here, we inserted two of them in the same function, `α`, to change the maximum value of the curve (at r = 1) and `cut`, the distance `r` where the function equals zero. In the image below, we can see the curve and how it changes with different values of `α` and `cut`:
 
 !!! INSERT FIGURE