keep working

app/checking.f90

@@ -0,0 +1,74 @@
+program main
+   use yaeos, only: pr, R
+   use yaeos, only: Groups, UNIFAC, setup_unifac
+   use yaeos, only: UNIQUAC, setup_uniquac
+
+   integer, parameter :: nc = 3, ng = 4
+
+   type(UNIFAC) :: unif
+   type(UNIQUAC) :: uniq
+   type(Groups) :: molecules(nc)
+
+   real(pr) :: Ge, Gen(nc), GeT, GeT2, GeTn(nc), Gen2(nc, nc)
+   real(pr) :: Ge_i, Gen_i(nc), GeT_i, GeT2_i, GeTn_i(nc), Gen2_i(nc, nc)
+   real(pr) :: ln_gammas(nc), b(nc, nc), rs(nc), qs(nc)
+
+   real(pr) :: n(nc), T, n_t
+
+   T = 150
+   n = [0.0_pr, 70.0_pr, 10.0_pr]
+   n_t = sum(n)
+
+   ! ! Ethane [CH3]
+   molecules(1)%groups_ids = [1]
+   molecules(1)%number_of_groups = [2]
+
+   ! ! Ethanol [CH3, CH2, OH]
+   molecules(2)%groups_ids = [1, 2, 14]
+   molecules(2)%number_of_groups = [1, 1, 1]
+
+   ! ! Methylamine [H3C-NH2]
+   molecules(3)%groups_ids = [28]
+   molecules(3)%number_of_groups = [1]
+
+   unif = setup_unifac(molecules)
+
+   ! ===========================================================================
+   ! UNIFAC
+   ! ---------------------------------------------------------------------------
+   call unif%excess_gibbs(n, T, Ge, GeT, GeT2, Gen, GeTn, Gen2)
+
+   print *, "UNIFAC:"
+   print *, "Ge: ", Ge
+   print *, "GeT: ", GeT
+   print *, "GeT2: ", GeT2
+   print *, "Gen: ", Gen
+   print *, "GeTn: ", GeTn
+   print *, "Gen2: ", Gen2
+
+   ! ===========================================================================
+   ! UNIQUAC
+   ! ---------------------------------------------------------------------------
+   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]
+
+   uniq = setup_uniquac(qs, rs, bij=b)
+
+   call uniq%excess_gibbs(n, T, Ge, GeT, GeT2, Gen, GeTn, Gen2)
+
+   print *, "UNIQUAC:"
+   print *, "Ge: ", Ge
+   print *, "GeT: ", GeT
+   print *, "GeT2: ", GeT2
+   print *, "Gen: ", Gen
+   print *, "GeTn: ", GeTn
+   print *, "Gen2: ", Gen2
+
+end program main

doc/page/usage/excessmodels/unifaclv.md

@@ -2,6 +2,8 @@
 title: UNIFAC-LV
 ---
 
+[TOC]
+
 # UNIFAC
 
 [[UNIFAC]] (UNIQUAC Functional-group Activity Coefficients) is an Excess Gibbs
@@ -15,25 +17,21 @@ $$ \frac{G^E}{RT} = \frac{G^{E,r}}{RT} + \frac{G^{E,c}}{RT} $$
 
 Being:
 
-- Combinatorial: Accounts for the size and shape of the molecules. The involved
-equations can be checked in the API documentation:
-[[yaeos__models_ge_group_contribution_unifac:Ge_combinatorial]]
+- Combinatorial: Accounts for the size and shape of the molecules.
 
-- Residual: Accounts for the energy interactions between different functional
-groups. The involved equations can be checked in the API documentation:
-[[yaeos__models_ge_group_contribution_unifac:Ge_residual]]
+- Residual: Accounts for the energy interactions between different functional groups.
 
 Each substance of a mixture modeled with [[UNIFAC]] must be represented by a
 list a functional groups and other list with the ocurrence of each functional
 group on the substance. The list of the functional groups culd be accesed on
 the DDBST web page:
 [https://www.ddbst.com/published-parameters-unifac.html](https://www.ddbst.com/published-parameters-unifac.html)
 
-# Examples
-## Calculating activity coefficients
-We can instantiate a [[UNIFAC]] model with a mixture  ethanol-water and
-evaluate the logarithm of activity coefficients of the model for a 0.5 mole
-fraction of each, and a temperature of 298.0 K.
+## Examples
+### Calculating activity coefficients
+We can instantiate a [[UNIFAC]] model with a mixture ethanol-water and evaluate
+the logarithm of activity coefficients of the model for a 0.5 mole fraction of
+each, and a temperature of 298.0 K.
 
 ```fortran
 use yaeos__constants, only: pr

doc/page/usage/excessmodels/uniquac.md

@@ -29,6 +29,24 @@ $$
 e_{lk}{T^2}
 $$
 
+Some of the model's terms can be simplified to reduce the complexity of the
+derivatives. Also, allows the model to be evaluated in mole vector `n` where
+some of the composition are equal to zero.
+
+$$\frac{\phi_k}{x_k} = \frac{n_T r_k}{\sum_l r_l n_l}$$
+
+$$\frac{\theta_k}{\phi_k} = \frac{q_k \sum_l r_l n_l}{r_k \sum_l q_l n_l}$$
+
+Being \(n_T \) the total number of moles in the system. The expression for the
+Excess Gibbs free energy can be rewritten as:
+
+$$ 
+\frac{G^E}{RT} = \sum_k n_k \ln \left(\frac{n_T r_k}{\sum_l r_l n_l} \right)
++ \frac{z}{2}\sum_k q_k n_k \ln \left(\frac{q_k \sum_l r_l n_l}{r_k \sum_l q_l
+n_l} \right) - \sum_k q_k n_k \ln\left(\sum_l \theta_l \tau_{lk} \right) 
+$$
+
+
 ## Temperature derivatives
 
 \(\qquad \tau_{lk}:\)
@@ -138,6 +156,19 @@ $$
 {\tau}_{lk}}
 $$
 
+Nueva:
+$$
+\frac{\partial \frac{G^E}{RT}}{\partial n_i} = \ln 
+\left(\frac{n_T r_i}{\sum_l r_l n_l} \right) + \sum_k n_k 
+\left(\frac{r_k \sum_l r_l n_l - n_T r_k r_i}{(\sum_l r_l n_l)^2}\right) +
+\frac{z}{2} q_i n_i \ln \left(\frac{q_i \sum_l r_l n_l}{r_i \sum_l q_l n_l} 
+\right) - {q}_{i}
+\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}}
+$$
+
+
 
 Differentiating each term of the first compositional derivative respect to
 \(n_j\) we get:

src/models/excess_gibbs/uniquac.f90

@@ -330,8 +330,8 @@ subroutine excess_gibbs(self, n, T, Ge, GeT, GeT2, Gen, GeTn, Gen2)
       ! ------------------------------------------------------------------------
       ! Combinatorial term
       Ge_comb = ( &
-         sum(n * log(phik / xk)) &
-         + self%z / 2.0_pr * sum(n * self%qs * log(thetak / phik)) &
+         sum(n * log(n_tot * self%rs / sum_nr)) &
+         + self%z / 2.0_pr * sum(n * self%qs * log(self%qs * sum_nr / self%rs / sum_nq)) &
          )
 
       ! Residual term
@@ -360,9 +360,10 @@ subroutine excess_gibbs(self, n, T, Ge, GeT, GeT2, Gen, GeTn, Gen2)
          do i=1,nc
             ! Combinatorial term
             Gen_comb(i) = (&
-               log(phik(i) / xk(i)) + sum(n * (dphik_dn(:,i) / phik - dxk_dni(:,i) / xk)) &
-               + self%z / 2.0_pr * self%qs(i) * log(thetak(i) / phik(i)) &
-               + self%z / 2.0_pr * sum(n * self%qs * (dthetak_dni(:,i) / thetak - dphik_dn(:,i) / phik)) &
+               log(n_tot * self%rs(i) / sum_nr) &
+               + sum(n * (self%rs * sum_nr - n_tot * self%rs(i)* self%rs) / sum_nr**2) &
+               + self%z / 2.0_pr * self%qs(i) * n(i) * log(self%qs(i) * sum_nr / self%rs(i) / sum_nq) &
+               + self%z / 2.0_pr * sum(n * self%qs * (self%rs(i) * sum_nq - self%qs(i) * sum_nr) / sum_nq / sum_nr) &
                )
 
             ! Residual term