Customize equilibrium function (with spatial derivatives)?

Thanks for your help in advance,
I am trying to implement the Cahn-Hilliard equation in PALABOS to simulate
multiphase flow with high density ratios (see Banari et. al implementation for
reference “An improved two-phase Lattice Boltzmann model for high density
ratios: application to wave breaking”). The equilibrium function is expressed as
(equation 21):


The first and final terms are not problematic, but the second term has a second
order spatial derivative with respect to $\phi$ (equation 2 in Banari et al.):


From my understanding all dynamics objects have a computeEquilibrium
method that operates cell-wise and always takes these arguments (from

virtual void computeEquilibrium(plint iPop, T rhoBar,
Array<T,Descriptor<T>::d> const& j, T jSqr, T thetaBar=T())

My approach was to have an additional term in this expression (e.g. T lap_rho)
that passes the laplacian as an argument to computeEquilibrium, but I am
unable to find where the computeEquilibrium method is actually called.
Put another way, I am trying to modify the call to computeEquilibirium so that
it can have the spatial derivatives as another argument.
I assume this call is integrated into the collideAndStream method of the
MultiBlockLattice, but extending the MultiBlockLattice seems to be going against
the recommendations of the User Guide.

Thanks again for help.