use yaeos, only: Groups, setup_unifac, UNIFAC
use yaeos__consistency_gemodel, only: numeric_ge_derivatives

type(UNIFAC) :: model

integer, parameter :: nc = 4, ng = 4

type(Groups) :: molecules(nc)

real(pr) :: Ge, Gen(nc), GeT, GeT2, GeTn(nc), Gen2(nc, nc)
real(pr) :: Ge_n, Gen_n(nc), GeT_n, GeT2_n, GeTn_n(nc), Gen2_n(nc, nc)
real(pr) :: ln_gammas(nc)

real(pr) :: n(nc), T
real(pr) :: dt, dn

T = 303.15
n = [400.0, 100.0, 300.0, 200.0] ! always test with sum(n) > 1

dt = 0.1_pr
dn = 0.1_pr

! Hexane [CH3, CH2]
molecules(1)%groups_ids = [1, 2]
molecules(1)%number_of_groups = [2, 4]

! Ethanol [CH3, CH2, OH]
molecules(2)%groups_ids = [1, 2, 14]
molecules(2)%number_of_groups = [1, 1, 1]

! Toluene [ACH, ACCH3]
molecules(3)%groups_ids = [9, 11]
molecules(3)%number_of_groups = [5, 1]

! Cyclohexane [CH2]
molecules(4)%groups_ids = [2]
molecules(4)%number_of_groups = [6]

model = setup_unifac(molecules)

! =====================================================================
! Call analytic derivatives
! ---------------------------------------------------------------------
call model%excess_gibbs(n, T, Ge, GeT, GeT2, Gen, GeTn, Gen2)

! =====================================================================
! Call numeric derivatives
! ---------------------------------------------------------------------
call numeric_ge_derivatives(model, n, T, dn, dt, Ge=Ge_n, GeT=GeT_n)
call numeric_ge_derivatives(model, n, T, dn, dt, Ge=Ge_n, Gen=Gen_n)
call numeric_ge_derivatives(model, n, T, dn, dt, Ge=Ge_n, GeT2=GeT2_n)
call numeric_ge_derivatives(model, n, T, dn, dt, Ge=Ge_n, GeTn=GeTn_n)
call numeric_ge_derivatives(model, n, T, dn, dt, Ge=Ge_n, Gen2=Gen2_n)