UNIQUAC done

c_interface/yaeos_c.f90

@@ -29,6 +29,7 @@ module yaeos_c
  ! GeMoels
  public :: nrtl
  public :: unifac_vle
+ public :: uniquac
  public :: ln_gamma
 
  ! Thermoprops
@@ -93,6 +94,28 @@ subroutine unifac_vle(id, nc, ngs, g_ids, g_v)
  call extend_ge_models_list(id)
  end subroutine
 
+ subroutine uniquac(id, qs, rs, aij, bij, cij, dij, eij)
+ use yaeos, only: setup_uniquac
+ integer(c_int), intent(out) :: id
+ real(c_double), intent(in) :: qs(:)
+ !! Molecule's relative areas \(Q_i\)
+ real(c_double), intent(in) :: rs(size(qs))
+ !! Molecule's relative volumes \(R_i\)
+ real(c_double), intent(in) :: aij(size(qs),size(qs))
+ !! Interaction parameters matrix \(a_{ij}\)
+ real(c_double), intent(in) :: bij(size(qs),size(qs))
+ !! Interaction parameters matrix \(b_{ij}\)
+ real(c_double), intent(in) :: cij(size(qs),size(qs))
+ !! Interaction parameters matrix \(c_{ij}\)
+ real(c_double), intent(in) :: dij(size(qs),size(qs))
+ !! Interaction parameters matrix \(d_{ij}\)
+ real(c_double), intent(in) :: eij(size(qs),size(qs))
+ !! Interaction parameters matrix \(e_{ij}\)
+
+ ge_model = setup_uniquac(qs, rs, aij, bij, cij, dij, eij)
+ call extend_ge_models_list(id)
+ end subroutine
+
  subroutine extend_ge_models_list(id)
  !! Find the first available model container and allocate the model
  !! there. Then return the found id.

doc/page/usage/excessmodels/uniquac.md

@@ -11,7 +11,7 @@ UNIQUAC (**uni**versal **qua**si**c**hemical) Excess Gibbs free energy model.
 $$ 
 \frac{G^E}{RT} = \sum_k n_k \ln\frac{\phi_k}{x_k}
 + \frac{z}{2}\sum_k q_k n_k \ln\frac{\theta_k}{\phi_k}
-- \sum_k q_k n_k \ln\left(\sum_l \theta_l \tau_{kl} \right)
+- \sum_k q_k n_k \ln\left(\sum_l \theta_l \tau_{lk} \right)
 $$
 
 With:
@@ -48,29 +48,29 @@ $$
 
 $$
 \frac{\partial G^E}{\partial T} = \frac{G^E}{T} - RT \sum_k q_k n_k \frac{
-\sum_l \theta_l \frac{\partial \tau_{kl}}{\partial T}}{\sum_l \theta_l
-\tau_{kl}}
+\sum_l \theta_l \frac{\partial \tau_{lk}}{\partial T}}{\sum_l \theta_l
+\tau_{lk}}
 $$
 
 $$
 \frac{\partial G^E}{\partial T^2} = -R\left[T \sum_k q_k n_k
-\left(\frac{(\sum_l \theta_l \frac{\partial^2 \tau_{kl}}{\partial T^2})}{\sum_l
-\theta_l \tau_{kl}}
-- \frac{(\sum_l \theta_l \frac{\partial \tau_{kl}}{\partial T})^2}{(\sum_l
-\theta_l \tau_{kl})^2}\right) + 2\left(\sum_k q_k n_k \frac{(\sum_l \theta_l
-\frac{\partial \tau_{kl}}{\partial T} )}{\sum_l \theta_l \tau_{kl}}\right)
+\left(\frac{(\sum_l \theta_l \frac{\partial^2 \tau_{lk}}{\partial T^2})}{\sum_l
+\theta_l \tau_{lk}}
+- \frac{(\sum_l \theta_l \frac{\partial \tau_{lk}}{\partial T})^2}{(\sum_l
+\theta_l \tau_{lk})^2}\right) + 2\left(\sum_k q_k n_k \frac{(\sum_l \theta_l
+\frac{\partial \tau_{lk}}{\partial T} )}{\sum_l \theta_l \tau_{lk}}\right)
 \right]
 $$
 
 ## Cross temperature-compositional derivative
 
 $$
 \frac{\partial^2 G^E}{\partial n_i \partial T} = \frac{1}{T} \frac{\partial
-G^E}{\partial n_i} - R T \left(q_i \frac{\sum_l \theta_l \frac{d \tau_{il}}{d
-T}}{\sum_l \theta_l \tau_{il}} + \sum_k q_k n_k \left(\frac{(\sum_l \frac{d
-\theta_l}{d n_i} \frac{d \tau_{kl}}{d T})(\sum_l \theta_l \tau_{kl}) - (\sum_l
-\theta_l \frac{d \tau_{kl}}{d T})(\sum_l \frac{d \theta_l}{d n_i}
-\tau_{kl})}{(\sum_l \theta_l \tau_{kl})^2} \right) \right)
+G^E}{\partial n_i} - R T \left(q_i \frac{\sum_l \theta_l \frac{d \tau_{li}}{d
+T}}{\sum_l \theta_l \tau_{li}} + \sum_k q_k n_k \left(\frac{(\sum_l \frac{d
+\theta_l}{d n_i} \frac{d \tau_{lk}}{d T})(\sum_l \theta_l \tau_{lk}) - (\sum_l
+\theta_l \frac{d \tau_{lk}}{d T})(\sum_l \frac{d \theta_l}{d n_i}
+\tau_{lk})}{(\sum_l \theta_l \tau_{lk})^2} \right) \right)
 $$
 
 
@@ -133,7 +133,7 @@ $$
 \frac{z}{2}{q}_{i}\ln{\left(\frac{\theta_{i}}{\phi_{i}} \right)} + \frac{z}{2}
 \sum_k {n}_{k} {q}_{k} \left(\frac{\frac{d \theta_{k}}{d {n}_{i}}}{\theta_k} -
 \frac{\frac{d \phi_{k}}{d {n}_{i}}}{\phi_k} \right) - {q}_{i}
-\ln{\left(\sum_l \theta_{l} {\tau}_{il} \right)} - \sum_k {n}_{k} {q}_{k}
+\ln{\left(\sum_l \theta_{l} {\tau}_{li} \right)} - \sum_k {n}_{k} {q}_{k}
 \frac{\sum_l \frac{d \theta_{l}}{d {n}_{i}} {\tau}_{lk}}{\sum_l \theta_{l}
 {\tau}_{lk}}
 $$
@@ -172,17 +172,17 @@ $$
 $$
 
 $$
-- q_i \left( \frac{\sum_l \frac{d \theta_l}{d n_j} \tau_{il}}{\sum_l \theta_l 
-\tau_{il}} \right)
+- q_i \left( \frac{\sum_l \frac{d \theta_l}{d n_j} \tau_{li}}{\sum_l \theta_l 
+\tau_{li}} \right)
 $$
 
 $$
-- {q}_{j} \frac{\sum_l \frac{d \theta_{l}}{d {n}_{i}} {\tau}_{jl}}{\sum_l
-\theta_{l}{\tau}_{jl}} - \sum_k {n}_{k} {q}_{k} \frac{\left(\sum_l
-\frac{d^2\theta_l}{dn_idn_j} \tau_{kl} \right) \left(\sum_l
-\theta_l \tau_{kl} \right) - \left(\sum_l \frac{d\theta_l}{dn_i}
-\tau_{kl} \right) \left(\sum_l \frac{d\theta_l}{dn_j} \tau_{kl}
-\right)}{(\sum_l \theta_{l} {\tau}_{kl})^2}
+- {q}_{j} \frac{\sum_l \frac{d \theta_{l}}{d {n}_{i}} {\tau}_{lj}}{\sum_l
+\theta_{l}{\tau}_{lj}} - \sum_k {n}_{k} {q}_{k} \frac{\left(\sum_l
+\frac{d^2\theta_l}{dn_idn_j} \tau_{lk} \right) \left(\sum_l
+\theta_l \tau_{lk} \right) - \left(\sum_l \frac{d\theta_l}{dn_i}
+\tau_{lk} \right) \left(\sum_l \frac{d\theta_l}{dn_j} \tau_{lk}
+\right)}{(\sum_l \theta_{l} {\tau}_{lk})^2}
 $$
 
 ## Examples
@@ -217,14 +217,13 @@ qs = [1.4_pr, 1.972_pr, 2.4_pr]
 T = 298.15_pr
 
 ! Calculate bij from DUij. We need -DU/R to get bij
-b = reshape(&
-[0.0_pr, -526.02_pr, -309.64_pr, &
-318.06_pr, 0.0_pr, 91.532_pr, &
--1325.1_pr, -302.57_pr, 0.0_pr], [nc, nc])
+b(1,:) = [0.0_pr, -526.02_pr, -309.64_pr]
+b(2,:) = [318.06_pr, 0.0_pr, 91.532_pr]
+b(3,:) = [-1325.1_pr, -302.57_pr, 0.0_pr]
 
 model = setup_uniquac(qs, rs, bij=b)
 
-n = [0.8_pr, 0.1_pr, 0.2_pr]
+n = [2.0_pr, 2.0_pr, 8.0_pr]
 
 call model%ln_activity_coefficient(n, T, ln_gammas)
 

python/yaeos/models/excess_gibbs/uniquac.py

@@ -55,17 +55,28 @@ class UNIQUAC(GeModel):
  -------
  .. code-block:: python
 
- from yaeos import UNIFACVLE
+ import numpy as np
 
- # Groups for water and ethanol
- water = {16: 1}
- ethanol = {1: 1, 2: 1, 14: 1}
+ from yaeos import UNIQUAC
 
- groups = [water, ethanol]
+ b = np.array(
+ [
+ [0.0, -526.02, -309.64],
+ [318.06, 0.0, 91.532],
+ [-1325.1, -302.57, 0.0],
+ ]
+ )
+
+ rs = np.array([0.92, 2.1055, 3.1878])
+ qs = np.array([1.4, 1.972, 2.4])
+
+ t = 298.15
+
+ model = UNIQUAC(qs, rs, bij=b)
 
- model = UNIFAVLE(groups)
+ n = np.array([2.0, 2.0, 8.0])
 
- model.ln_gamma([0.5, 0.5], 298.15)
+ gammas = np.exp(model.ln_gamma(n, t)) # [8.856, 0.860, 1.425]
  """
 
  def __init__(

src/models/excess_gibbs/uniquac.f90

@@ -46,6 +46,41 @@ module yaeos__models_ge_uniquac
  !! d_{lk}T + e_{lk}{T^2}
  !! \]
  !!
+ !! # Example
+ !!
+ !! ```fortran
+ !! use yaeos, only: pr, setup_uniquac, UNIQUAC
+ !!
+ !! integer, parameter :: nc = 3
+ !!
+ !! real(pr) :: rs(nc), qs(nc)
+ !! real(pr) :: b(nc, nc)
+ !! real(pr) :: n(nc)
+ !!
+ !! real(pr) :: ln_gammas(nc), T
+ !!
+ !! type(UNIQUAC) :: model
+ !!
+ !! rs = [0.92_pr, 2.1055_pr, 3.1878_pr]
+ !! qs = [1.4_pr, 1.972_pr, 2.4_pr]
+ !!
+ !! T = 298.15_pr
+ !!
+ !! ! Calculate bij from DUij. We need -DU/R to get bij
+ !! b(1, :) = [0.0_pr, -526.02_pr, -309.64_pr]
+ !! b(2, :) = [318.06_pr, 0.0_pr, 91.532_pr]
+ !! b(3, :) = [-1325.1_pr, -302.57_pr, 0.0_pr]
+ !!
+ !! model = setup_uniquac(qs, rs, bij=b)
+ !!
+ !! n = [0.8_pr, 0.1_pr, 0.2_pr]
+ !!
+ !! call model%ln_activity_coefficient(n, T, ln_gammas)
+ !!
+ !! print *, exp(ln_gammas) ! [8.856, 0.860, 1.425]
+ !!
+ !! ```
+ !!
  real(pr) :: z = 10.0_pr
  !! Model coordination number
  real(pr), allocatable :: qs(:)
@@ -303,7 +338,7 @@ subroutine excess_gibbs(self, n, T, Ge, GeT, GeT2, Gen, GeTn, Gen2)
  sum_thetal_tau_lk = 0.0_pr
 
  do k=1,nc
- sum_thetal_tau_lk(k) = sum(thetak * tau(k,:))
+ sum_thetal_tau_lk(k) = sum(thetak * tau(:,k))
  end do
 
  Ge_res = -sum(n * self%qs * log(sum_thetal_tau_lk))
@@ -315,7 +350,7 @@ subroutine excess_gibbs(self, n, T, Ge, GeT, GeT2, Gen, GeTn, Gen2)
  ! ------------------------------------------------------------------------
  if (dn .or. dtn .or. dn2) then
  do concurrent (k=1:nc, i=1:nc)
- sum_dtheta_l_tau_lk(i,k) = sum(dthetak_dni(:,i) * tau(k,:))
+ sum_dtheta_l_tau_lk(i,k) = sum(dthetak_dni(:,i) * tau(:,k))
  end do
  end if
 
@@ -363,14 +398,14 @@ subroutine excess_gibbs(self, n, T, Ge, GeT, GeT2, Gen, GeTn, Gen2)
  ) &
  )
 
- trm5 = -self%qs(i) * (sum(dthetak_dni(:,j) * tau(i,:)) / sum_thetal_tau_lk(i))
+ trm5 = -self%qs(i) * (sum(dthetak_dni(:,j) * tau(:,i)) / sum_thetal_tau_lk(i))
 
  do k=1,nc
- sum_d2theta_tau_lk(k) = sum(d2thetak_dnidnj(:,i,j) * tau(k,:))
+ sum_d2theta_tau_lk(k) = sum(d2thetak_dnidnj(:,i,j) * tau(:,k))
  end do
 
  trm6 = (&
- -self%qs(j) * (sum(dthetak_dni(:,i) * tau(j,:)) / sum_thetal_tau_lk(j)) &
+ -self%qs(j) * (sum(dthetak_dni(:,i) * tau(:,j)) / sum_thetal_tau_lk(j)) &
  -sum(self%qs * n * (&
  sum_d2theta_tau_lk * sum_thetal_tau_lk &
  - sum_dtheta_l_tau_lk(i,:) * sum_dtheta_l_tau_lk(j,:) &
@@ -386,7 +421,7 @@ subroutine excess_gibbs(self, n, T, Ge, GeT, GeT2, Gen, GeTn, Gen2)
  sum_theta_l_dtau_lk = 0.0_pr
 
  do k=1,nc
- sum_theta_l_dtau_lk(k) = sum(thetak * dtau(k,:))
+ sum_theta_l_dtau_lk(k) = sum(thetak * dtau(:,k))
  end do
  end if
 
@@ -400,7 +435,7 @@ subroutine excess_gibbs(self, n, T, Ge, GeT, GeT2, Gen, GeTn, Gen2)
  sum_theta_l_d2tau_lk = 0.0_pr
 
  do k=1,nc
- sum_theta_l_d2tau_lk(k) = sum(thetak * d2tau(k,:))
+ sum_theta_l_d2tau_lk(k) = sum(thetak * d2tau(:,k))
  end do
 
  diff_aux = (&
@@ -417,7 +452,7 @@ subroutine excess_gibbs(self, n, T, Ge, GeT, GeT2, Gen, GeTn, Gen2)
  ! Cross derivative Tn
  if (dtn) then
  do concurrent (k=1:nc, i=1:nc)
- sum_dtheta_l_dtau_lk(i,k) = sum(dthetak_dni(:,i) * dtau(k,:))
+ sum_dtheta_l_dtau_lk(i,k) = sum(dthetak_dni(:,i) * dtau(:,k))
  end do
  end if
 
@@ -498,41 +533,7 @@ end subroutine taus
  type(UNIQUAC) function setup_uniquac(qs, rs, aij, bij, cij, dij, eij)
  !! Instantiate a UNIQUAC model.
  !!
- !! # Example
- !!
- !! ```fortran
- !! use yaeos, only: pr, setup_uniquac, UNIQUAC
- !!
- !! integer, parameter :: nc = 3
- !!
- !! real(pr) :: rs(nc), qs(nc)
- !! real(pr) :: b(nc, nc)
- !! real(pr) :: n(nc)
- !!
- !! real(pr) :: ln_gammas(nc), T
- !!
- !! type(UNIQUAC) :: model
- !!
- !! rs = [0.92_pr, 2.1055_pr, 3.1878_pr]
- !! qs = [1.4_pr, 1.972_pr, 2.4_pr]
- !!
- !! T = 298.15_pr
- !!
- !! ! Calculate bij from DUij. We need -DU/R to get bij
- !! b = reshape(&
- !! [0.0_pr, -526.02_pr, -309.64_pr, &
- !! 318.06_pr, 0.0_pr, 91.532_pr, &
- !! -1325.1_pr, -302.57_pr, 0.0_pr], [nc, nc])
- !!
- !! model = setup_uniquac(qs, rs, bij=b)
- !!
- !! n = [0.8_pr, 0.1_pr, 0.2_pr]
- !!
- !! call model%ln_activity_coefficient(n, T, ln_gammas)
- !!
- !! print *, exp(ln_gammas) ! [8.856, 0.860, 1.425]
- !!
- !! ```
+ !! Non provided interaction parameters are set to zero matrices.
  !!
  real(pr), intent(in) :: qs(:)
  !! Molecule's relative volumes \(Q_i\)