Reweight Function massive rework

README.md

@@ -4,6 +4,11 @@
 [![codecov](https://codecov.io/gh/m3g/MolSimToolkit.jl/branch/main/graph/badge.svg)](https://codecov.io/gh/m3g/MolSimToolkit.jl)
 [![Aqua QA](https://JuliaTesting.github.io/Aqua.jl/dev/assets/badge.svg)](https://github.com/JuliaTesting/Aqua.jl)
 
+# MyDevMolSimToolkit
+
+[MolSimToolkit.jl](https://github.com/m3g/MolSimToolkit.jl) provides a set of tools to 
+analyse molecular dynamics simulations, and a framework for the development of custom
+analysis tools. This repository is a experimental version of the package for Work in Progress functionalities
 # MolSimToolkit.jl
 
 **MolSimToolkit.jl** is a Julia package for analyzing molecular dynamics (MD) simulations. It provides fast, flexible tools for trajectory analysis, structural alignment, secondary structure mapping, and more. The toolkit is designed to be lightweight and extensible, making it easy to build custom analysis workflows.

docs/src/Reweighting.md

@@ -39,40 +39,46 @@ 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
@@ -125,6 +131,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]

src/MolSimToolkit.jl

@@ -25,6 +25,7 @@ using Reexport: @reexport
 using ProgressMeter: Progress, next!, @showprogress
 using Statistics: mean
 using Printf: @sprintf
+using OrderedCollections
 
 export wrap, wrap_to_first
 export distances