## Definitions

Within `Clapeyron`

, it is possible to create your own model/equation of state whilst still using all of the property estimation tools we provide. In order to do this, you must create the model. We provide a macro to make it easy to set it up. If you wish to create a new model called CustomEoS, you just need call the `@newmodel`

macro with three parameters:

`@newmodel CustomEoS CustomEosModel CustomEosParam`

- Model struct name
- Model abstract type
- Model parameters struct

We will now give a brief overview of each of these aspects.

### Model struct name

This is the concrete type, which is simply the name of the equation of state which will be used to generate the model. Within `Clapeyron`

, we've tried to keep these names as obvious as possible (*e.g.* `vdW`

, `PR`

, `PCSAFT`

, `SAFTVRMie`

, `GERG2008`

, *etc.*). For your own model, this can be whatever you want as long is the identification makes sense to you.

### Model abstract type

In Julia, objects can only be sub-types of abstract types (and not concrete types), which will allow you to inherit the behaviour of the parent(s). In order to maintain a hierarchy of models, we highly encourage you to create an abstract type for your model that is a subtype one of these families of EoS, and dispatch on this newly created abstract type instead of on your model's concrete type. This is not enforced, but we suffix all abstract types in this category with `Model`

as a convention. The abstract type that is at the top of the hierarchy is `EoSModel`

. From this parent, we branch into more-specific EoS sub-types:

`SAFTModel`

: These are the models which have three parameters in common: segment size, $\sigma$, potential depth, $\epsilon$, and number of segments, $m$. All other SAFT-type models branch from this parent (`PCSAFTModel`

,`SAFTVRMieModel`

,`softSAFTModel`

,*etc.*).`CubicModel`

: These are the models whose parameters can be obtained from the critical temperature and pressure. With the exception of`CPAModel`

, all cubics have a common structure where one can re-arrange the equation for the pressure as a third-order polynomial. As such, we define a subtype of`CubicModel`

,`ABCubicModel`

(*e.g.*`vdWModel`

,`RKModel`

,`SRKModel`

,`PRModel`

).`EmpiricHelmholtzModel`

: These are the high-accuracy, multi-parameter models for specific species or systems (*e.g.*`GERG2008Model`

,`IAPWS95Model`

). There is no general structure to the models and they are treated as self-contained.`IdealModel`

: Often overlooked, these models supplement the`SAFTModel`

and`CubicModel`

by providing the ideal contribution. Whilst the parameters and structure aren't usually the same between ideal models, this is unnecessary as the equation for the pressure is always $pV=Nk_\mathrm{B}T$ .

For example, if we wish to create a new EoS model called `CustomEoS`

, we will create an abstract type `CustomEoSModel`

, that is a sub-type of another abstract type, say `CubicModel`

(but it could also inherit from something further down the hierarchy, like `vdWModel`

):

`abstract type CustomEoSModel <: CubicModel end`

For models that are sub-types of `SAFTModel`

or `CubicModel`

, most methods will be instantly compatible because methods will be able to make use of a set of the EoS parameters for finding the initial guesses. If your model is not of either of these types, it will be necessary to define a few additional functions:

`lb_volume(model::CustomEoSModel,T,z)`

: This must output the smallest possible value of the volume for your particular model using the model parameters. In SAFT equations, this is equivalent to a packing fraction of one ($\eta=1$) whilst in cubics, it is equivalent to the`b`

parameter.`T_scale(model::CustomEoSModel,z)`

: This must output the temperature scaling for your model using the model parameters. In SAFT equations, this is usually the potential depth whilst in cubics it is the critical temperature.`p_scale(model::CustomEoSModel,z)`

: This must output the pressure scaling for your model using the model parameters. In cubics, this is the critical pressure whilst in SAFT we use a more complicated definition using the segment sizes and potential depths.

### Model parameters struct

The parameters for a particular system are all stored within a struct that is a subtype of `EoSParam`

. By convention, we suffix these with `Param`

. These structs should contain the model parameters, which comprise objects of types `SingleParam{T}`

, `PairParam{T}`

, and `AssocParam{T}`

, where `T`

is usually a base type (`Float64`

, `Integer`

, `String`

, etc). Below is an example of a generic param struct for a SAFT and cubic model.

```
struct GenericSAFTParam <: EoSParam
Mw::SingleParam{Float64}
segment::SingleParam{Float64}
sigma::PairParam{Float64}
epsilon::PairParam{Float64}
epsilon_assoc::AssocParam{Float64}
bondvol::AssocParam{Float64}
end
```

```
struct GenericCubicParam <: EoSParam
Tc::SingleParam{Float64}
pc::SingleParam{Float64}
Mw::SingleParam{Float64}
a::PairParam{Float64}
b::PairParam{Float64}
end
```

The parameters are wrapped in one of the following structs:

`SingleParam{T}`

: These are the parameters associated with a pure species or indexed by a single index ($i$). For example, the molecular weight, number of segments, critical temperature and pressure.`PairParam{T}`

: These are the parameters associated with a pair of species or indexed by two indices ($i$ and $j$). If the two indices are the same ($i=j$), they refer to the like species, otherwise ($i\neq j$), they refer to the unlike species. For example, the cubic`a`

and`b`

parameters, the SAFT segment size (`sigma`

) and potential depth (`epsilon`

) parameters.`AssocParam{T}`

: There are the parameters associated with both a pair of species and association sites (see background documentation for what these are). They are indexed by four indices (the species $i$ and $j$, and the sites $a$ and $b$). For example, potential depth of the association interaction (`epsilon_assoc`

) and the bonding volume (`bondvol`

).

If there exists a model parameter struct that contains exactly the same parameters as the one that you are about to create, you can also directly use that existing struct.

## PC-SAFT Example

Once all the above has been defined, we are ready to build our own model. Let us imagine we are trying to implement `PCSAFT`

.

- We first define the Model name, type and parameters:

```
# Defining an abstract type for this model type
abstract type PCSAFTModel <: SAFTModel end
# Defining the parameters used by the model
struct PCSAFTParam <: EoSParam
Mw::SingleParam{Float64}
segment::SingleParam{Float64}
sigma::PairParam{Float64}
epsilon::PairParam{Float64}
epsilon_assoc::AssocParam{Float64}
bondvol::AssocParam{Float64}
end
# Creating a model struct called PCSAFT, which is a sub-type of PCSAFTModel, and uses parameters defined in PCSAFTParam
@newmodel PCSAFT PCSAFTModel PCSAFTParam
```

- We define some functions to inform how the model obtains their parameters, what are their references, and if we need to perform any reescaling, combining or renaming of input parameters:

```
# Use this function if you want to provide default parameters for your new model.
# Those locations are relative to the default database location.
Clapeyron.default_locations(::Type{PCSAFT}) = ["PCSAFT_parameters/PCSAFT_like.csv","PCSAFT_parameters/PCSAFT_unlike.csv","PCSAFT_parameters/PCSAFT_assoc.csv"]
# Although optional, it's generally good practise to cite your models!
# We use DOI references in general.
Clapeyron.default_references(::Type{PCSAFT}) = ["10.1021/ie0003887", "10.1021/ie010954d"]
# Use this function if you need to perform any transformations, from the input parameters, to the final ones:
function Clapeyron.transform_params(::Type{PCSAFT},params)
# params is a Dict{String,ClapeyronParam}
sigma = params["sigma"]
# We reescale the sigma values
sigma.values .*= 1e-10
# If k is not provided, it will be not be considered
k = get(params,"k",nothing)
# In some cases, we may not have the unlike parameters and will need to use combining rules. You can also define your own combining rules for this.
sigma = sigma_LorentzBerthelot(sigma)
epsilon = epsilon_LorentzBerthelot(params["epsilon"], k)
# Now, we put the parameters in the params dict
params["sigma"] = sigma
params["epsilon"] = epsilon
# And, finally, we return the dictionary with the final values
return params
end
```

The `@newmodel`

macro, in conjunction with the functions defined above, gets lowered to the following julia code:

```
struct PCSAFT{T<:IdealModel} <: PCSAFTModel
components::Array{String,1} # Components vector. the model gets it's length from here
sites::SiteParam # Parameter struct containing the sites and their amounts
params::PCSAFTParam # Struct specified in the macro
idealmodel::T # Model for the ideal part
assoc_options::Clapeyron.AssocOptions # Options for the calculation of the association helmholtz contribution
references::Array{String,1} # DOI references
end
# Constructor generated by the macro:
function PCSAFT(components; # Vector{String} containing a list of components
idealmodel=BasicIdeal, # Default ideal model
userlocations=String[], # User locations for the main model
ideal_userlocations=String[], # User locations for the ideal model
verbose=false, # Verbosity
assoc_options = AssocOptions()) # Options for the calculation of the association contribution
# Obtain a Dict of parameters. We pass in custom locations through the optional parameter userlocations.
params = getparams(components,default_locations(model); userlocations=userlocations, verbose=verbose)
# We perform the specified transformation of parameters.
transform_params(PCSAFT,params)
# We perform mixing of the association parameters
# depending on the options contained in assoc_options
sigma = params["sigma"]
epsilon_assoc = params["epsilon_assoc"]
bondvol = params["bondvol"]
# Combining rules for association. if you want to perform cross-association mixing, check the AssocOptions docs
bondvol,epsilon_assoc = assoc_mix(bondvol,epsilon_assoc,sigma,assoc_options)
# Now we can create the parameter struct that we have defined.
packagedparams = PCSAFTParam(params["Mw"], params["segment"], sigma, params["epsilon"], epsilon_assoc, bondvol)
references = default_references(PCSAFT)
#we initialize the ideal model provided.
init_idealmodel = Clapeyron.init_model(idealmodel,components,ideal_userlocations,verbose)
# Build the model.
model = PCSAFT(components,sites,packagedparams,init_idealmodel,assoc_options,references)
# Return the PCSAFT object that you have just created.
return model
end
```

You can, of course, not use the macro, if your model depends itself on other models (like `CPA`

depending on a cubic model),or requires more arguments than the ones generated by the macro (like mixing rules for cubics that require an activity model)

Define all of the model equations. We encourage you to use the full range of Unicode characters where it makes your code clearer to read!

As convention, the first four arguments should be

`model`

,`V`

,`T`

and`z`

; any other variables should come after.If we obey that convention, we may use the

`@f`

macro, which automatically substitutes the first four parameters for compactness. For example,`@f(func,i,j)`

is equivalent to calling`func(model,V,T,z,i,j)`

.Clapeyron obtains all the properties of a model by differenciating the total helmoltz energy (

`eos`

) or the residual helmoltz energy (`eos_res`

).`eos`

and`eos_res`

themselves are defined in terms of the reduced ideal helmholtz energy (`a_res`

). In this case, we are going to define`a_res`

for our own model:`function Clapeyron.a_res(model::PCSAFTModel, V, T, z) return @f(a_hc) + @f(a_disp) + @f(a_assoc) end function a_hc(model::PCSAFTModel, V, T, z) x = z/∑(z) m = model.params.segment.values m̄ = ∑(x .* m) return m̄*@f(a_hs) - ∑(x[i]*(m[i]-1)*log(@f(g_hs,i,i)) for i ∈ @comps) end function d(model::PCSAFTModel, V, T, z, i) ϵii = model.params.epsilon.values[i,i] σii = model.params.sigma.values[i,i] return σii * (1 - 0.12exp(-3ϵii/T)) end function ζ(model::PCSAFTModel, V, T, z, n) ∑z = ∑(z) x = z * (one(∑z)/∑z) m = model.params.segment.values res = N_A*∑z*π/6/V * ∑((x[i]*m[i]*@f(d,i)^n for i ∈ @comps)) end function g_hs(model::PCSAFTModel, V, T, z, i, j) di = @f(d,i) dj = @f(d,j) ζ2 = @f(ζ,2) ζ3 = @f(ζ,3) return 1/(1-ζ3) + di*dj/(di+dj)*3ζ2/(1-ζ3)^2 + (di*dj/(di+dj))^2*2ζ2^2/(1-ζ3)^3 end function a_hs(model::PCSAFTModel, V, T, z) ζ0 = @f(ζ,0) ζ1 = @f(ζ,1) ζ2 = @f(ζ,2) ζ3 = @f(ζ,3) return 1/ζ0 * (3ζ1*ζ2/(1-ζ3) + ζ2^3/(ζ3*(1-ζ3)^2) + (ζ2^3/ζ3^2-ζ0)*log(1-ζ3)) end # INSERT REST OF CODE`

With all the above defined in a single script, we can save the file as

`PCSAFT.jl`

and then include it in our jupyter notebooks (for example) and use the model with all of our existing method:`include("PCSAFT.jl") model = PCSAFT(["carbon dioxide"]) p = 20e6 T = range(290,460,length=200) Cp = isobaric_heat_capacity.(model, p, T) (T_c, p_c, V_c) = crit_pure(model) T_sat = range(220,T_c,length=200) (p_sat, V_l_sat, V_v_sat) = saturation_pressure(model,T_sat)`

## sPC-SAFT Example

Instead of developing an entirely new model, some of us may want to modify or extend an existing one. `sPCSAFT`

is an example where we want to modify parts of regular `PCSAFT`

but keep the rest the same. We can do this in a very succinct way making this new model a sub-type of the abstract type associated with another model.

When we define the model type,

`sPCSAFT`

is a sub-type of`PCSAFT`

:`abstract type sPCSAFTModel <: PCSAFTModel end`

Since the parameters are the same, we can just use the same model params when creating the model:

`@newmodel sPCSAFT sPCSAFTModel PCSAFTParam`

This may not be the case if we're extending a model (

*e.g.*if we're adding polar or ionic terms, we may need to define a new parameter struct to include the new parameters).When defining the model equations, we only need to write those that have been changed in

`sPCSAFT`

:`function a_hc(model::sPCSAFTModel, V, T, z) x = z/sum(z) m = model.params.segment.values m̄ = ∑(x .* m) return m̄*@f(a_hs) - (m̄-1)*log(@f(g_hs)) end function g_hs(model::sPCSAFTModel, V, T, z) η = @f(ζ,3) return (1-η/2)/(1-η)^3 end function a_hs(model::sPCSAFTModel, V, T, z) η = @f(ζ,3) return (4η-3η^2)/(1-η)^2 end`

The rest works exactly as it would with the `PCSAFT`

example.