ion pair simulation

Dear all,
I am calculating Internal energy using a derivation of an equation of state, and I need to calculate a temperature dependent parameter( which is total energy in the ideal gas state(F(T)):

Ein=∫Pin/rho^2 d rho +F(T)= RT{ (erho^2)+(frho)+(g*rho^4)}+F(T)

That F(T) is the ideal contribution of internal energy since it is equal to total energy at zero density
Also, the resulted Cp from the derivation of this Ein is:

Cp=R{ (arho^2)+(brho)+(c*rho^4)}+F’(T)
This Cp has F’(T) which is temperature derivation of F(T).
. I simulated bulk for the non-ideal part with 300 molecules of [BMIM][PF6], and it is totally ok. But for the ideal part(F(T) or F’(T)) at gas state, energies are not consistent and so meaningful, because as I calculate F(T) from experimental data, I found out that total energy should be negative, which is now positive in Lammps(for one ion pair), and smaller than what it is now. I used a single pair of the ion at the nonperiodic boundary condition, but the energies( bond, angle, dihedral, …) are not the same as bulk which I think they should be since they are intramolecular quantities(they are slightly bigger than what they should be as regard to article that I use)(also non-bonded intramolecular energies like VdW and electrostatic are slightly different). Also keeping the temperature constant is so difficult, but it is not a problem, and I can derive any F(T) or F’(T) with fitting them into a polynomial. I already used 4RT for the Cp as F’(T)(or ideal contribution of Cp), but it did not work, so I want to use Lammps to calculate energies at each temperature and use their energy as F(T), and finally, I would use the derivation of F(T) instead of F’(T)for the ideal contribution of Cp.

any suggestion or comment would be appreciated.

best regards