Arni/extract henry coeff

REQUIRE

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

src/Misc.jl

@@ -136,3 +136,81 @@ function read_atomic_masses()
 
     return atomic_masses
 end
+
+
+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
+        M0 = n[end] * 1.1
+        if isapprox(n[1], 0.0)
+            K0 = n[2]/p[2]
+        else
+            K0 = n[1]/p[1]
+        end
+        return Dict("M0" => M0, "K0" => K0)
+    elseif model == :henry
+        if isapprox(n[1], 0.0)
+            K0 = n[2]/p[2]
+        else
+            K0 = n[1]/p[1]
+        end
+        return Dict("K0" => K0)
+    else
+        error("Model not available. Currently only `:langmuir` and `:henry` are available")
+    end
+end
+
+
+
+"""
+    params = fit_isotherm(df, pressure_col_name, loading_col_name, model)
+
+Takes in a DataFrame `df` containing an isotherm. Will try to fit a model to the data.
+Available models are `:henry` and `:langmuir`
+The Langmuir model is in 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 partial pressure of the gas.
+`params` will contain a fit to M and K
+The Henry model takes the following form:
+N = KP
+and `params` will contain K
+
+# 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 the isotherm
+
+# 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 "K".
+"""
+function fit_isotherm(df::DataFrame, pressure_col_name::Symbol, loading_col_name::Symbol, model::Symbol; henry_tol::Float64 = 0.25)
+    sort!(df, [pressure_col_name])
+    n = df[loading_col_name]
+    p = df[pressure_col_name]
+    if ! isapprox(p[1], 0.0, atol=0.01)
+        prepend!(n, 0.0)
+        prepend!(p, 0.0)
+    end
+    θ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
+        min_mse = Inf
+        objective_function_henry(θ) = return sum([(n[i] - θ[1] * p[i])^2 for i = 1:length(n)])
+        res = optimize(objective_function_henry, [θ0["K0"]] * 1.1, LBFGS())
+        K = res.minimizer
+        mse = res.minimum / length(n)
+        if mse < min_mse
+            min_mse = mse
+        end
+        return Dict("K" => K[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_isotherm,
 
     # Crystal.jl
     Framework, read_crystal_structure_file, remove_overlapping_atoms_and_charges,

test/isotherm.csv

@@ -0,0 +1,9 @@
+Pressure(bar),Adsorption(mmol/g)
+0.0,0.0
+0.560,0.534
+1.333,1.295
+2.370,2.204
+3.551,3.153
+4.830,4.080
+6.128,4.951
+7.244,5.616

test/misc_test.jl

@@ -0,0 +1,29 @@
+module Misc_Test
+
+using PorousMaterials
+using Test
+using DataFrames
+using CSV
+using Optim
+
+@testset "Misc Tests" begin
+    df = CSV.read("isotherm.csv")
+    p_col_name, l_col_name = names(df)
+    henry = fit_isotherm(df[1:3,:], p_col_name, l_col_name, :henry)["K"]
+    @test isapprox(henry, 0.9688044280548711)
+
+    P = range(0, stop=1, length=100)
+    M = 23.0
+    K = 11.9
+    N = (M*K.*P)./(1 .+ K.*P)
+    N = [N[1], N[11], N[21], N[31], N[41], N[51], N[61], N[71], N[81]]
+    P = [P[1], P[11], P[21], P[31], P[41], P[51], P[61], P[71], P[81]]
+    new_df = DataFrame(Pressure = P, Loading = N)
+    x = fit_isotherm(new_df, Symbol("Pressure"), Symbol("Loading"), :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")