Theoretical vs Simulated coexistence curve

I am trying to replicate some results from Lycett-Brown (2015) (see ref. at the end of this post), specifically, Figure 8 with the Guo forcing scheme. According to this work (Eq. 12), by using a pseudopotential function psi = exp(-1/rho), Guo forcing scheme coexistence curve matches exactly the Maxwell construction rule for the planar interface problem. I have implemented this by modifying the shanChenMultiPhase example. Basically, I modified the original data processor to compute Fx and Fy and saves these values to the cell external field related to force. As for the dynamics, I am using GuoExternalForceBGKDynamics to include the effects of the force computed by the data processor. I was able to replicate the interparticle potential described by using PsiShanChen94 with psi0 = rho0 = 1. I am initializing the liquid and gas densities with the expected theorical values, and using equation (37) of Kharmiani et al (2019). I am running 50000 time steps, and results seem to achieve a steady state way before this. For the time being, I am using a 5 x 200 grid. I was expecting densities to be very close to Maxwell construction for this case, however simulation results achieve equilibrium with values way off from the expected ones. I have already revised the data processor (implementation, level values, etc…), dynamics collision model (equilibrium and Guo forcing equations inside Palabos code), and tried different strategies such as turning off the default roundoff policy (it seemed to me that default roundoff makes sense when densities are very close to 1. However for this case where densities varies from 0.2 to 4.3, I thought this could be causing some unexpected effects), distinct forcing schemes, and distinct interparticle potential. However, results are always off from the expected ones (see figure). I am not sure what I am missing.
I would be glad if someone could help me to sort this out. I am sharing my setup if someone would like to reproduce what I am trying to simulate. Also, I have attached the modified process method. I guess it may be some silly mistake, but I have spend some time double checking everything, and I was not able to identify what is causing this behavior.

Simulation setup:
Dynamics: GuoExternalForceBGKdynamics (omega = 1)
Grid: 5 x 200
Boundary conditions: periodic in all directions
Data processor (level 1): modified ShanChenDataProcessor (force saved to external field)
Interparticle Potential: ShanChen94 (psi0 =1, rho0 = 1)
Lattice descriptor: ForcedD2Q9Descriptor
Initialization: Equation (37) in Kharmiani et al (2019). Basically, a smooth tanh function from rho_g and rho_l, with rho_l at central nodes of simulation.
Initialization values: expected rho_g and rho_l from Maxwell curve

G;rho_g;rho_l
-13.0;0.21673836314073247;4.29731824006538
-12.198436585561522;0.2482758470927612;3.897263279087027
-11.396873171123042;0.2848090756365337;3.4952278102450767
-10.595309756684564;0.3284049756844069;3.0891423216904688
-9.793746342246086;0.38294339821251316;2.6751869215784394
-8.992182927807608;0.45666264519052735;2.2452529087031614
-8.190619513369128;0.5725512710446172;1.776556833741509

template<typename T, template<typename U> class Descriptor>
void SCForceProcessor2D<T,Descriptor>::process (
        Box2D domain, BlockLattice2D<T,Descriptor>& lattice )
{
    // Short-hand notation for the lattice descriptor
    typedef Descriptor<T> D;
    // Handle to external scalars
    enum {
        forceOffset  = D::ExternalField::forceBeginsAt
    };

    plint nx = domain.getNx() + 2;  // Include a one-cell boundary
    plint ny = domain.getNy() + 2;  // Include a one-cell boundary
    plint offsetX = domain.x0-1;
    plint offsetY = domain.y0-1;
    ScalarField2D<T> psiField(nx,ny);
    
    for (plint iX=domain.x0-1; iX<=domain.x1+1; ++iX) {
        for (plint iY=domain.y0-1; iY<=domain.y1+1; ++iY) {
            // Get "intelligent" value of density through cell object, to account
            //   for the fact that the density value can be user-defined, for example
            //   on boundaries.
            Cell<T,Descriptor>& cell = lattice.get(iX,iY);
            T rho = cell.computeDensity();
            // Evaluate potential function psi.
            psiField.get(iX-offsetX, iY-offsetY) = Psi->compute(rho);
        }
    }

    // Compute the interparticle forces
    for (plint iX=domain.x0; iX<=domain.x1; ++iX) {
        for (plint iY=domain.y0; iY<=domain.y1; ++iY) {
            Array<T,D::d> rhoContribution;
            rhoContribution.resetToZero();
            // Compute the term \sum_i ( t_i psi(x+c_i,t) c_i )
            for (plint iPop = 0; iPop < D::q; ++iPop) {
                plint nextX = iX + D::c[iPop][0];
                plint nextY = iY + D::c[iPop][1];
                T psi = psiField.get(nextX-offsetX, nextY-offsetY);
                for (int iD = 0; iD < D::d; ++iD) {
                    rhoContribution[iD] += D::t[iPop] * psi * D::c[iPop][iD];
                }
            }

            Cell<T,Descriptor>& cell = lattice.get(iX,iY);
            T *force = cell.getExternal(forceOffset);
            for (int iD = 0; iD < D::d; ++iD) {
                T psi = psiField.get(iX-offsetX, iY-offsetY);
                force[iD] = -G * psi * rhoContribution[iD];
            }
        }

References:
Lycett-Brown, Daniel, and Kai H. Luo. “Improved forcing scheme in pseudopotential lattice Boltzmann methods for multiphase flow at arbitrarily high density ratios.” Physical Review E 91, no. 2 (2015): 023305.

Kharmiani, Soroush Fallah, Hamid Niazmand, and Mohammad Passandideh-Fard. “An Alternative High-Density Ratio Pseudo-potential Lattice Boltzmann Model with Surface Tension Adjustment Capability.” Journal of Statistical Physics 175, no. 1 (2019): 47-70.