The pressure drop with periodic boundary condition can be calculated as following. Suppose you need to achieve u_final with your simulation. I mean that your final velocity which is involved in Poiseuille velocity distribution for the flow inside plates equals to

u(x)= grad Pressure / visc (x^2/2-H^2/8), where H is Ny in our case - the width of channel. That means that in the center the velocity is gradPressure/viscH^2/8. If it’s equal to u_final we can obtain formula that gradPressure =rho08visc*ufinal/(Ny^2)

Now the force representation is involved. If you’d like to represent the force in the simplest case, like Fc_i/c_s^2, then to get force in terms of pressure you need to multiply this on c_s^2. After that you can obtain that gradPressure=Fc_i.

Therefore force_i term which we need to add to the right hand side of the LBEquation will be gradPressure/6, where 6 is the number of population in which gradient of pressure directs or flow goes (D2Q9 has 6 velocities with x projection).

Finally fpois equals to rho0* 8.0/6.0 * visc * uf / dfloat(ny) / dfloat(ny) and you add it to the ith population as fpois*c_i/c_s^2. You will obtain the u_final as the center velocity in the channel.

Hopefully my math is OK and it can help you a bit,
Alex