add tols func

src/helper_funcs.jl

@@ -24,6 +24,15 @@ function nearest_mirror(r_ij, box_sizes)
     return [r_x,r_y,r_z] 
 end
 
+function apply_tols!(arr, tol)
+    Threads.@threads for i in eachindex(arr)
+        if abs(arr[i]) < tol
+            arr[i] = 0.0
+        end
+    end
+    return arr
+end
+
 # """
 # nᵗʰ triangular number
 # """

src/mcc.jl

@@ -5,15 +5,17 @@ export mcc3
 Converts third order forces constants, `Ψ` into third order modal coupling constants (MCC).
 Does not divide MCC calculation into smaller chunks. This might exhaust GPU memory.
 """
-function mcc3(Ψ::CuArray{Float32, 3}, phi::CuArray{Float32, 2})
+function mcc3(Ψ::CuArray{Float32, 3}, phi::CuArray{Float32, 2}, tol)
 
     K3 = CUDA.zeros(Float32, size(Ψ));
 
     @tensor begin
         K3[n,m,l] = Ψ[i,j,k] * phi[i,n] * phi[j,m] * phi[k,l]
     end
 
-    return Array(K3)
+    K3_CPU = Array(K3)
+    return apply_tols!(K3_CPU, tol)
+
 end
 
 """
@@ -23,7 +25,7 @@ parameter `block_size` specifies problem size when calculating the MCC. For exam
 The `phi` matrix will be automatically truncated to adjust for this.
     Try to maximize `block_size` and make it a power of 2.
 """
-function mcc3(Ψ::CuArray{Float32, 3}, phi::CuArray{Float32, 2}, block_size::Int)
+function mcc3(Ψ::CuArray{Float32, 3}, phi::CuArray{Float32, 2}, block_size::Int, tol)
 
     @assert size(phi)[1] % block_size == 0
     @assert block_size > 0
@@ -53,7 +55,7 @@ function mcc3(Ψ::CuArray{Float32, 3}, phi::CuArray{Float32, 2}, block_size::Int
     #TODO this calculates more than it needs to atm, can only calculate unique part and rotate 
     #TODO to enforce symmetry, kinda annoying to implement, will matter more when system > GPU RAM
 
-    return K3_CPU
+    return apply_tols!(K3_CPU, tol)
 end
 
 

src/third_order.jl

@@ -69,19 +69,14 @@ function third_order_IFC(sys::SuperCellSystem{D}, pot::PairPotential, tol) where
     end
 
     #Apply tolerances
-    for i in eachindex(Ψ)
-        if abs(Ψ[i]) < tol
-            Ψ[i] = 0.0
-        end
-    end
+    Ψ = apply_tols!(Ψ, tol)
 
     #Give proper units
     Ψ_unit = unit(pot.ϵ / pot.σ^3)
 
     return ThirdOrderMatrix(Ψ, Ψ_unit, tol)
 end
 
-
 function mass_weight_third_order!(Ψ::ThirdOrderMatrix, masses::AbstractVector)
     N_modes = size(Ψ)[1]
     D = Int(N_modes/length(masses))