This notebook computes the chemical equilibrium composition of the Haber process (N₂ + 3H₂ ⇌ 2NH₃) using a non-stoichiometric Gibbs energy minimization approach, following the method of Binous and Bellagi (Engineering Reports, 2022). Temperature-dependent molar enthalpy and entropy for hydrogen, nitrogen, and ammonia are computed from NIST WebBook Shomate equation parameters, from which the species’ partial molar Gibbs free energies, GjH2, GjN2, and GjNH3, are assembled as functions of temperature. Each species’ Gibbs free energy, G[T, P, j], further incorporates a pressure correction term via R T Log[P/P1], referencing a standard pressure of 1 bar.
The total Gibbs free energy of the reacting mixture, Fitness[{n1, n2, n3}], is built by summing each species’ contribution — its chemical potential plus an ideal-mixing entropy term, R Tn Log[n[j]/Sum[n[k]]] — weighted by its mole number. Since the minimization is performed over unconstrained mole numbers rather than a single reaction-extent variable, atomic mass balance (fixed total nitrogen and hydrogen atoms, based on a 3:1 H₂:N₂ feed) is enforced indirectly through two quadratic penalty terms with coefficient 10⁹, which strongly discourage any candidate solution from violating the balance equations while keeping the objective differentiable and search-friendly for a swarm-based optimizer.
A custom Particle Swarm Optimization routine, psoOptimize, searches this penalized objective over a 3-dimensional mole-number space bounded by FitnessLB and FitnessUB. Using 1000 particles over 300 iterations with linearly decaying inertia weight (0.9 → 0.4) and standard cognitive/social coefficients (c1 = c2 = 1.42), the swarm converges from an initial loss near 487 (dominated by constraint violation) down to a stable minimum of approximately −88.61 by iteration ~120, after which the loss plateaus, indicating convergence to a global optimum candidate at n1 ≈ 0.363, n2 ≈ 0.121, n3 ≈ 0.258.
The PSO solution is refined and validated against FindMinimum, a local Newton-based solver seeded at the PSO’s best-found point, which converges to a more precise optimum of −88.6433 at n1 = 0.410118, n2 = 0.136706, n3 = 0.226588, with residual mass-balance errors on the order of 10⁻⁸ — confirming the constraints are satisfied to high precision.
To visualize convergence behavior meaningfully despite the loss trajectory crossing from positive to negative values, the plot tracks the gap between each iteration’s loss and the final best loss (convergence - gbestscore) on a log scale, revealing the characteristic PSO convergence pattern: a steep initial descent through several orders of magnitude followed by a long, flat tail as the swarm fine-tunes near the optimum.
The total Gibbs free energy of the reacting mixture, Fitness[{n1, n2, n3}], is built by summing each species’ contribution — its chemical potential plus an ideal-mixing entropy term, R Tn Log[n[j]/Sum[n[k]]] — weighted by its mole number. Since the minimization is performed over unconstrained mole numbers rather than a single reaction-extent variable, atomic mass balance (fixed total nitrogen and hydrogen atoms, based on a 3:1 H₂:N₂ feed) is enforced indirectly through two quadratic penalty terms with coefficient 10⁹, which strongly discourage any candidate solution from violating the balance equations while keeping the objective differentiable and search-friendly for a swarm-based optimizer.
A custom Particle Swarm Optimization routine, psoOptimize, searches this penalized objective over a 3-dimensional mole-number space bounded by FitnessLB and FitnessUB. Using 1000 particles over 300 iterations with linearly decaying inertia weight (0.9 → 0.4) and standard cognitive/social coefficients (c1 = c2 = 1.42), the swarm converges from an initial loss near 487 (dominated by constraint violation) down to a stable minimum of approximately −88.61 by iteration ~120, after which the loss plateaus, indicating convergence to a global optimum candidate at n1 ≈ 0.363, n2 ≈ 0.121, n3 ≈ 0.258.
The PSO solution is refined and validated against FindMinimum, a local Newton-based solver seeded at the PSO’s best-found point, which converges to a more precise optimum of −88.6433 at n1 = 0.410118, n2 = 0.136706, n3 = 0.226588, with residual mass-balance errors on the order of 10⁻⁸ — confirming the constraints are satisfied to high precision.
To visualize convergence behavior meaningfully despite the loss trajectory crossing from positive to negative values, the plot tracks the gap between each iteration’s loss and the final best loss (convergence - gbestscore) on a log scale, revealing the characteristic PSO convergence pattern: a steep initial descent through several orders of magnitude followed by a long, flat tail as the swarm fine-tunes near the optimum.
Thermodynamic data obtained from the NIST WebBook : https://webbook.nist.gov/
Hydrogen Section
Hydrogen Section
WBH2={33.066178,-11.363417,11.432816,-2.772874,-0.158558,-9.980797,172.707974,0.0};(*ShomateEquationParametersforH2*)t=Tn/1000;THH2={t,t^2/2,t^3/3,t^4/4,-1/t,1,0,-1};TSH2={Log[t],t,t^2/2,t^3/3,-1/(2*t^2),0,1,0};HH2=Sum[(WBH2*THH2)[[i]],{i,1,8}];SH2=Sum[(WBH2*TSH2/1000)[[i]],{i,1,8}];H298H2=0.000;GjH2[Tn_]:=H298H2+HH2-Tn*SH2;
Nitrogen Section
Nitrogen Section
WBN2={19.50583,19.88705,-8.598535,1.369784,0.527601,-4.935202,212.3900,0.0};(*ShomateEquationParametersforN2*)t=Tn/1000;THN2={t,t^2/2,t^3/3,t^4/4,-1/t,1,0,-1};TSN2={Log[t],t,t^2/2,t^3/3,-1/(2*t^2),0,1,0};HN2=Sum[(WBN2*THN2)[[i]],{i,1,8}];SN2=Sum[(WBN2*TSN2/1000)[[i]],{i,1,8}];H298N2=0.000;GjN2[Tn_]:=H298N2+HN2-Tn*SN2;
Ammonia Section
Ammonia Section
WBNH3={19.99563,49.77119,-15.37599,1.921168,0.189174,-53.30667,203.8591,-45.89806};(*ShomateEquationParametersforNH3*)t=Tn/1000;THNH3={t,t^2/2,t^3/3,t^4/4,-1/t,1,0,-1};TSNH3={Log[t],t,t^2/2,t^3/3,-1/(2*t^2),0,1,0};HNH3=Sum[(WBNH3*THNH3)[[i]],{i,1,8}];SNH3=Sum[(WBNH3*TSNH3/1000)[[i]],{i,1,8}];H298NH3=-45.940;GjNH3[Tn_]:=H298NH3+HNH3-Tn*SNH3;
Components: H2 (1) N2 (2) NH3 (3)
G[x_,P_,1]:=GjH2[x]+RxLog[P/P1];G[x_,P_,2]:=GjN2[x]+RxLog[P/P1];G[x_,P_,3]:=GjNH3[x]+RxLog[P/P1];
In[]:=
P1=1;R=8.31410^-3;
Mole fractions
In[]:=
y[i_]:=n[i]/Sum[n[j],{j,1,3}];
Total Gibbs free energy assuming contribution of Residual Gibbs free energy is negligible.
Fitness function is sum of total Gibbs free energy and penalty functions (to Enforce Exact Mass Balance Constraints)
Fitness function is sum of total Gibbs free energy and penalty functions (to Enforce Exact Mass Balance Constraints)
Fitness[{n1_,n2_,n3_}]:=Sum[n[j]*(G[Tn,P,j]+RTnLog[n[j]/Sum[n[k],{k,1,3}]]),{j,1,3}]+10^9(2n1+3n3-2×0.75)^2+10^9(2n2+n3-2×0.25)^2/.n[1]->n1/.n[2]->n2/.n[3]->n3;
3D problem
SwarmDim=3;
Upper and lower bounds for the search domain
FitnessLB={0.05,0.05,0.05};FitnessUB={0.5,0.2,0.3};
PSO parameters
c1=1.42;c2=1.42;w=0.72984;
In[]:=
Clear[n1,n2,n3]
Without loss of generality : pressure and temperature are taken equal to 600 bar and 800 K
P=600;Tn=800;
In[]:=
Fitness[{n1,n2,n3}]
Out[]=
1000000000+1000000000+n1-70.3897+6.6512Log+n2-119.22+6.6512Log+n3-169.121+6.6512Log
2
(-0.5+2n2+n3)
2
(-1.5+2n1+3n3)
n1
n1+n2+n3
n2
n1+n2+n3
n3
n1+n2+n3
PSO code
Clip each column of a matrix against its own scalar bound
clipMatrix[mat_,lbv_,ubv_]:=Transpose[MapThread[Clip[#1,{#2,#3}]&,{Transpose[mat],lbv,ubv}]];clampMatrix[mat_,vmaxv_]:=Transpose[MapThread[Clip[#1,{-#2,#2}]&,{Transpose[mat],vmaxv}]];
PSO Optimizer
psoOptimize[lossFunc_,dim_,lb_,ub_,nParticles_:30,nIter_:100,wMax_:0.9,wMin_:0.4,c1_:1.42,c2_:1.42,vClampFrac_:0.2]:=Module[{particles,velocities,pbest,pbestScore,gbest,gbestScore,convergence={},r1,r2,scores,gbestIdx,vmax,wIter},vmax=vClampFrac*(ub-lb);particles=Transpose@Table[RandomReal[{lb[[i]],ub[[i]]},nParticles],{i,dim}];velocities=ConstantArray[0.,{nParticles,dim}];pbest=particles;pbestScore=lossFunc/@particles;gbestIdx=First@Ordering[pbestScore,1];gbestScore=pbestScore[[gbestIdx]];gbest=pbest[[gbestIdx]];Do[wIter=wMax-(wMax-wMin)*(j-1)/(nIter-1);r1=RandomReal[1,{nParticles,dim}];r2=RandomReal[1,{nParticles,dim}];velocities=wIter*velocities+c1*r1*(pbest-particles)+c2*r2*(ConstantArray[gbest,nParticles]-particles);velocities=clampMatrix[velocities,vmax];particles+=velocities;particles=clipMatrix[particles,lb,ub];scores=lossFunc/@particles;pbest=MapThread[If[#3<#2,#1,#4]&,{particles,pbestScore,scores,pbest}];pbestScore=MapThread[Min,{pbestScore,scores}];gbestIdx=First@Ordering[pbestScore,1];If[pbestScore[[gbestIdx]]<gbestScore,gbestScore=pbestScore[[gbestIdx]];gbest=pbest[[gbestIdx]];];AppendTo[convergence,gbestScore];If[Mod[j,10]==0,Print["PSO Iter ",j," Loss = ",gbestScore]];,{j,nIter}];{gbest,gbestScore,convergence}];
In[]:=
SeedRandom[42];
PSO run
result=psoOptimize[Fitness,SwarmDim,FitnessLB,FitnessUB,1000,300];{gbest,gbestScore,convergence}=result;
PSO Iter 10 Loss = 487.153
PSO Iter 20 Loss = 487.153
PSO Iter 30 Loss = 110.034
PSO Iter 40 Loss = -12.6574
PSO Iter 50 Loss = -42.2514
PSO Iter 60 Loss = -79.917
PSO Iter 70 Loss = -88.1893
PSO Iter 80 Loss = -88.1893
PSO Iter 90 Loss = -88.6024
PSO Iter 100 Loss = -88.6024
PSO Iter 110 Loss = -88.6077
PSO Iter 120 Loss = -88.6085
PSO Iter 130 Loss = -88.6085
PSO Iter 140 Loss = -88.6085
PSO Iter 150 Loss = -88.6085
PSO Iter 160 Loss = -88.6085
PSO Iter 170 Loss = -88.6085
PSO Iter 180 Loss = -88.6085
PSO Iter 190 Loss = -88.6085
PSO Iter 200 Loss = -88.6085
PSO Iter 210 Loss = -88.6085
PSO Iter 220 Loss = -88.6085
PSO Iter 230 Loss = -88.6085
PSO Iter 240 Loss = -88.6085
PSO Iter 250 Loss = -88.6085
PSO Iter 260 Loss = -88.6085
PSO Iter 270 Loss = -88.6085
PSO Iter 280 Loss = -88.6085
PSO Iter 290 Loss = -88.6085
PSO Iter 300 Loss = -88.6085
In[]:=
gbest
Gap to the best-known solution at each iteration
Check material balance constraints
References
References
Housam Binous and Ahmed Bellagi, Calculations of complex chemical reaction equilibria using stoichiometric and non-stoichiometric approaches in combination with arc-length continuation, Engineering Reports, 4(10) , 2022.
CITE THIS NOTEBOOK
CITE THIS NOTEBOOK
Study of the Haber process using hybrid PSO (global search) + FindMinimum (local Newton-type refinement)
by Housam Binous
Wolfram Community, STAFF PICKS, July 3, 2026
https://community.wolfram.com/groups/-/m/t/3743669
by Housam Binous
Wolfram Community, STAFF PICKS, July 3, 2026
https://community.wolfram.com/groups/-/m/t/3743669