keep working tomorrow

app/cosito.f90

@@ -0,0 +1,48 @@
+program main
+   use yaeos, only: pr, R
+   use yaeos, only: Groups, setup_unifac, UNIFAC
+
+   type(UNIFAC) :: model
+
+   integer, parameter :: nc = 3, ng = 4
+
+   type(Groups) :: molecules(nc)
+
+   real(pr) :: n(nc), T, n_t
+
+   real(pr) :: He, HeT, Hen(nc)
+   real(pr) :: Se, SeT, Sen(nc)
+
+   T = 150
+   n = [20.0, 70.0, 10.0]
+   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]
+
+   ! setup UNIFAC model
+   model = setup_unifac(molecules)
+
+   call model%excess_enthalpy(n, T, He=He, HeT=HeT, Hen=Hen)
+   call model%excess_entropy(n, T, Se=Se, SeT=SeT, Sen=Sen)
+
+   print *, 'Excess enthalpy:', He / n_t, -8.12666342863807_pr
+   print *, 'Excess entropy:', Se / n_t, -0.03268447167877293_pr
+
+   print *, 'dHe/dT:', HeT / n_t, 0.05391608033744438_pr
+   print *, 'dSe/dT:', SeT / n_t, 0.0003594405355829625_pr
+
+   print *, 'dHe/dn:', Hen / n_t, [150.72393613,  -573.71657621, -4412.09526745]
+   print *, 'dSe/dn:', Sen / n_t, [-6.01538894_pr, -2.23972167_pr, -4.97564211_pr]
+
+
+end program main

coso.ipynb

@@ -0,0 +1,166 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": 8,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from thermo import UNIFAC\n",
+    "from thermo.unifac import UFIP, UFSG\n",
+    "import numpy as np\n",
+    "\n",
+    "\n",
+    "GE = UNIFAC.from_subgroups(\n",
+    "    chemgroups=[{1:2}, {1:1, 2:1, 14:1}, {28: 1}], \n",
+    "    T=150, \n",
+    "    xs=[0.2, 0.7, 0.1], \n",
+    "    version=0, \n",
+    "    interaction_data=UFIP, \n",
+    "    subgroups=UFSG\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 9,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "-0.03268447167877293"
+      ]
+     },
+     "execution_count": 9,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "GE.SE() / 100"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 10,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "0.0003594405355829625"
+      ]
+     },
+     "execution_count": 10,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "GE.dSE_dT() / 100"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 22,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "array([-6.01538894, -2.23972167, -4.97564211])"
+      ]
+     },
+     "execution_count": 22,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "np.array(GE.dSE_dxs())"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 14,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "-8.12666342863807"
+      ]
+     },
+     "execution_count": 14,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "GE.HE() / 100"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 15,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "0.05391608033744438"
+      ]
+     },
+     "execution_count": 15,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "GE.dHE_dT() / 100"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 25,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "array([  9.63390279,   2.38949767, -35.99428925])"
+      ]
+     },
+     "execution_count": 25,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "np.array(GE.dHE_dns()) / 100"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "yaeos",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.10.12"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}

src/models/excess_gibbs/ge_models.f90

@@ -9,21 +9,24 @@ module yaeos__models_ge
    contains
       procedure(excess_gibbs), deferred :: excess_gibbs
       procedure :: ln_activity_coefficient => ln_activity_coefficient
+      procedure :: excess_enthalpy => excess_enthalpy
+      procedure :: excess_entropy => excess_entropy
    end type
 
    abstract interface
-      subroutine excess_gibbs(self, n, t, Ge, GeT, GeT2, Gen, GeTn, Gen2)
-         !! Excess Gibbs and derivs procedure
+      subroutine excess_gibbs(self, n, T, Ge, GeT, GeT2, Gen, GeTn, Gen2)
+         !! Calculate Excess Gibbs and its derivatives.
+         !!
          import pr, GeModel
          class(GeModel), intent(in) :: self !! Model
          real(pr), intent(in) ::n(:) !! Moles vector
-         real(pr), intent(in) :: t !! Temperature [K]
-         real(pr), optional, intent(out) :: Ge !! Excess Gibbs
+         real(pr), intent(in) :: T !! Temperature [K]
+         real(pr), optional, intent(out) :: Ge !! Excess Gibbs free energy
          real(pr), optional, intent(out) :: GeT !! \(\frac{dG^E}{dT}\)
          real(pr), optional, intent(out) :: GeT2 !! \(\frac{d^2G^E}{dT^2}\)
-         real(pr), optional, intent(out) :: Gen(size(n))
-         real(pr), optional, intent(out) :: GeTn(size(n))
-         real(pr), optional, intent(out) :: Gen2(size(n), size(n))
+         real(pr), optional, intent(out) :: Gen(size(n)) !! \(\frac{dG}{dn}\)
+         real(pr), optional, intent(out) :: GeTn(size(n)) !! \(\frac{d^2G^E}{dTdn}\)
+         real(pr), optional, intent(out) :: Gen2(size(n), size(n)) !! \(\frac{d^2G^E}{dn^2}\)
       end subroutine
    end interface
 
@@ -40,4 +43,58 @@ subroutine ln_activity_coefficient(self, n, T, lngamma)
       call self%excess_gibbs(n, t, ge=ge, gen=dgedn)
       lngamma = dgedn/(R*T)
    end subroutine
+
+   subroutine excess_enthalpy(self, n, T, He, HeT, Hen)
+      !! Calculate Excess enthalpy and its derivatives.
+      !!
+      !! \[
+      !!  \frac{\partial \frac{G^E}{T}}{\partial T} = \frac{-H^E}{T^2}
+      !! \]
+      !!
+      !! \[
+      !!  H^E = G^E - T \frac{\partial G^E}{\partial T}
+      !! \]
+      !!
+      class(GeModel), intent(in) :: self !! Model
+      real(pr), intent(in) :: n(:) !! Moles vector
+      real(pr), intent(in) :: T !! Temperature [K]
+      real(pr), optional, intent(out) :: He !! Excess enthalpy
+      real(pr), optional, intent(out) :: HeT !! \(\frac{dH^E}{dT}\)
+      real(pr), optional, intent(out) :: Hen(:) !! \(\frac{dH^E}{dn}\)
+
+      real(pr) :: Ge, GeT, GeT2, Gen(size(n)), GeTn(size(n))
+
+      call self%excess_gibbs(&
+         n, T, Ge=Ge, GeT=GeT, GeT2=GeT2, Gen=Gen, GeTn=GeTn &
+         )
+
+      if (present(He)) He = Ge - T*GeT
+      if (present(HeT)) HeT = -T * GeT2
+      if (present(Hen)) Hen = Gen - T*GeTn
+   end subroutine excess_enthalpy
+
+   subroutine excess_entropy(self, n, T, Se, SeT, Sen)
+      !! Calculate Excess entropy and its derivatives.
+      !!
+      !! \[
+      !! S^E = \frac{H^E - G^E}{T}
+      !! \]
+      !!
+      class(GeModel), intent(in) :: self !! Model
+      real(pr), intent(in) :: n(:) !! Moles vector
+      real(pr), intent(in) :: T !! Temperature [K]
+      real(pr), optional, intent(out) :: Se !! Excess entropy
+      real(pr), optional, intent(out) :: SeT !! \(\frac{dS^E}{dT}\)
+      real(pr), optional, intent(out) :: Sen(:) !! \(\frac{dS^E}{dn}\)
+
+      real(pr) :: Ge, GeT, Gen(size(n))
+      real(pr) :: He, HeT, Hen(size(n))
+
+      call self%excess_gibbs(n, T, Ge=Ge, GeT=GeT, Gen=Gen)
+      call self%excess_enthalpy(n, T, He=He, HeT=HeT, Hen=Hen)
+
+      if (present(Se)) Se = (He - Ge) / T
+      if (present(SeT)) SeT = ((HeT - GeT) * T - He + Ge) / T**2
+      if (present(Sen)) Sen = (Hen - Gen) / T
+   end subroutine excess_entropy
 end module