Reweight Extension Modifications (Work in Progres)

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"

docs/make.jl

@@ -20,7 +20,7 @@ makedocs(
         "Developer zone" => "Developer.md",
         "Plotting style" => "plotting_style.md",
         "Experimental" => "Experimental.md",
-        " ∘ Simulation Reweighting" => "Reweighting.md",
+        " ∘ Simulation Reweighting" => "Reweighting.md"
     ],
 )
 deploydocs(

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/Reweighting/Reweight.jl

@@ -68,7 +68,7 @@
 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
 -------------
@@ -109,7 +109,7 @@
 function reweight(
     simulation::Simulation, 
     f_perturbation::Function, 
-    group_1::AbstractVector{<:Integer}; 
+    group_1::AbstractVector{<:Integer};
     cutoff::Real = 12.0, 
     k::Real = 1.0, 
     T::Real = 1.0
@@ -127,9 +127,9 @@
             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)
+        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
@@ -138,9 +138,9 @@
     simulation::Simulation, 
     f_perturbation::Function, 
     group_1::AbstractVector{<:Integer}, 
-    group_2::AbstractVector{<:Integer};     
+    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))
@@ -158,9 +158,9 @@
             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)
+        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

src/Reweighting/Reweighting.jl

@@ -2,6 +2,7 @@
 
 using CellListMap: ParticleSystem, map_pairwise, map_pairwise!
 using ..MolSimToolkit: Simulation, positions, unitcell
+using ..MolSimToolkit.MolecularMinimumDistances
 
 export reweight, lennard_jones_perturbation, poly_decay_perturbation, gaussian_decay_perturbation
 
@@ -23,7 +24,7 @@
 
     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
 
@@ -38,11 +39,17 @@
 
     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
@@ -62,9 +69,7 @@
 
     β = 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

src/Reweighting/perturbation_examples.jl

@@ -1,32 +1,32 @@
-#SETTING EQUATIONS FOR PERTUBATION
+#FUNCTIONS FOR PERTURBATIONS
 
 #Perturbing Lennard-Jones interactions
 function switching_func(r, switch, cut)
     switching = ((cut^2-r^2)^2*(cut^2-r^2-3*(switch^2-r^2)))/(cut^2-switch^2)^3
     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
+    return V_pert
 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 β")
+gaussian_decay_perturbation(r, alpha, beta) = beta > 0 ? alpha*exp(-beta*r^2) : error("you need to insert positive values for β")