Merge 59679f4 into adfae91

REQUIRE

@@ -10,3 +10,4 @@ ProgressMeter
 JLD2
 LightGraphs
 FileIO
+Optim

src/Misc.jl

@@ -136,3 +136,75 @@ function read_atomic_masses()
 
     return atomic_masses
 end
+
+# obtain a reasonable initial guess for the optimization route in fitting adsorption 
+#  models to adsorption isotherm data
+function _guess(df::DataFrame, pressure_col_name::Symbol, loading_col_name::Symbol, model::Symbol)
+    n = df[loading_col_name]
+    p = df[pressure_col_name]
+    if model == :langmuir || model == :henry
+        # Henry coefficient H guess as the secant line between origin and first data point 
+        #  of the adsorption isotherm
+        if isapprox(n[1], 0.0)
+            H0 = n[2] / p[2]
+        else
+            H0 = n[1] / p[1]
+        end
+        if model == :henry
+            return Dict("H0" => H0)
+        else
+            # saturation loading M is largest value of adsorption observed plus 10%
+            M0 = n[end] * 1.1
+            K0 = M0 * H0
+            return Dict("M0" => M0, "K0" => K0)
+        end
+    else
+        error("Model not available. Currently only `:langmuir` and `:henry` are available")
+    end
+end
+
+"""
+    params = fit_adsorption_isotherm(df, pressure_col_name, loading_col_name, model)
+
+Takes in a DataFrame `df` containing adsorption isotherm data and fits an analytical model 
+to the data to identify its parameters of best fit, returned as a dictionary.
+Available models are `:henry` and `:langmuir`
+
+The Henry model takes the following form:
+N = HP
+The identified Henry coefficient is `params["H"]`.
+
+The Langmuir model takes the following form:
+N = (MKP)/(1+KP)
+where N is the total adsorption, M is the maximum monolayer coverage, K is the Langmuir constant. and P is the pressure of the gas.
+
+# Arguments
+- `df::DataFrame`: The DataFrame containing the pressure and adsorption data for the isotherm
+- `pressure_col_name::Symbol`: The header of the pressure column. Can be found with `names(df)`
+- `loading_col_name::Symbol`: The header of the loading/adsorption column. Can be found with `names(df)`
+- `model::Symbol`: The model chosen to fit to the adsorption isotherm data
+
+# Returns
+- `params::Dict{AbstractString, Float64}`: A Dictionary with the parameters corresponding to each model along with the MSE of the fit. `:langmuir` contains "M" and "K". `:henry` contains "H".
+"""
+function fit_adsorption_isotherm(df::DataFrame, pressure_col_name::Symbol, 
+                                 loading_col_name::Symbol, model::Symbol)
+    sort!(df, [pressure_col_name])
+    n = df[loading_col_name]
+    p = df[pressure_col_name]
+    θ0 = _guess(df, pressure_col_name, loading_col_name, model)
+
+    if model == :langmuir
+        objective_function_langmuir(θ) = return sum([(n[i] - θ[1] * θ[2] * p[i] / (1 + θ[2] * p[i]))^2 for i = 1:length(n)])
+        res = optimize(objective_function_langmuir, [θ0["M0"], θ0["K0"]], LBFGS())
+        M, K = res.minimizer
+        mse = res.minimum / length(n)
+        return Dict("M" => M, "K" => K, "MSE" => mse)
+    elseif model == :henry
+        objective_function_henry(θ) = return sum([(n[i] - θ[1] * p[i])^2 for i = 1:length(n)])
+        res = optimize(objective_function_henry, [θ0["H0"]], LBFGS())
+        H = res.minimizer
+        mse = res.minimum / length(n)
+        return Dict("H" => H[1], "MSE" => mse)
+    end
+end

src/PorousMaterials.jl

@@ -14,6 +14,7 @@ using Printf
 using LinearAlgebra
 using LightGraphs
 using Distributed
+using Optim
 import Base.push!
 
 
@@ -106,7 +107,7 @@ export
     nearest_image!, nearest_r², nearest_r,
 
     # Misc.jl
-    read_xyz, read_cpk_colors, read_atomic_radii, write_xyz,
+    read_xyz, read_cpk_colors, read_atomic_radii, write_xyz, fit_adsorption_isotherm,
 
     # Crystal.jl
     Framework, read_crystal_structure_file, remove_overlapping_atoms_and_charges,

test/misc_test.jl

@@ -0,0 +1,32 @@
+module Misc_Test
+
+using PorousMaterials
+using Test
+using DataFrames
+using CSV
+using Optim
+
+@testset "Misc Tests" begin
+    ###
+    # adsorption isotherm tests
+    ###
+    # Henry
+    df = DataFrame(P = [0.0, 0.56, 1.333], N = [0.0, 0.534, 1.295])
+    henry = fit_adsorption_isotherm(df, :P, :N, :henry)["H"]
+    @test isapprox(henry, 0.9688044280548711)
+
+    # Langmuir
+    P = range(0, stop=1, length=100)
+    M = 23.0
+    K = 11.9
+    N = (M * K .* P) ./ (1 .+ K.*P)
+    ids_fit = [1, 11, 21, 31, 41, 51, 61, 71, 81]
+    df = DataFrame(P = P[ids_fit], N = N[ids_fit])
+    x = fit_adsorption_isotherm(df, :P, :N, :langmuir)
+    M2, K2 = x["M"], x["K"]
+    @test isapprox(M, M2)
+    @test isapprox(K, K2)
+end
+
+
+end

test/runtests.jl

@@ -16,4 +16,5 @@ include("gcmc_checkpoints_test.jl")
 include("henry_checkpoint_test.jl")
 include("grid_test.jl")
 include("path_test.jl")
+include("misc_test.jl")
 include("generic_rotation_test.jl")