Free energy simulation method is applied to calculate the virial coefficients of square-well (SW) fluids of variable well-width and square-well based dimer forming associating fluids. In this approach, Monte Carlo sampling is performed on a number of molecules equal to the order of integral, and configurations are weighted according to the absolute value of the integrand. An umbrella-sampling average yields the value of the cluster integral in reference to a known integral. By using this technique, we determine the virial coefficients up to B-6 for SW fluid with variable potential range from lambda = 1.25 to lambda = 3.0 and model associating fluids with different association strengths: epsilon(af) = 0.0, 8.0, 16.0 and 22.0. These calculated values for SW fluids are in good agreement with the literature. We examine these coefficients in the context of the virial equation of state (VEOS) of SW fluids. VEOS up to B-4 or up to B-6 describes the PVT behavior along the saturated vapor line better than the series that includes B-5. We used these coefficients to find the critical properties of SW fluids and compared with the literature values. Boyle temperature is also determined and is found to increase with the increase in the well-extent and associating strength. We also report Joule-Thomson inversion curve for Lennard-Jones fluid and SW fluids using different truncated VEOS and compared with that predicted from established EOS. (C) 2009 Elsevier B.V. All rights reserved.