The second, third, and fourth virial coefficients of hard chain molecules with number of segments up to 10 (up to 7 in the case of the fourth one) and the reduced distances L-*=1 and 0.5 were determined numerically. For fused hard sphere (FHS) molecules the enlarged fused hard sphere model is introduced which forms a connection to the hard convex geometry and makes it possible to determine the virial coefficients of FHS molecules from the expressions derived for hard convex bodies. Our numerical values of the virial coefficients together with data from literature are used to test the proposed method and to compare its results with values from Wertheim's theory, from its modified version and from the generalized Flory dimer approach. It is found that prediction of the second virial coefficient is very accurate (for L=0.5 the maximum deviation amounts 0.2 percent) and that our values of the third virial coefficient for higher members of the FHS family are superior to those from other considered methods. The model is successfully used to predict geometric characteristics of the hard models of n-butane conformers and to estimate their second virial coefficients. The equation of state of hard prolate spherocylinders in which the nonsphericity parameters of the enlarged FHS model are considered is used to predict the P-V-T behavior of hard dumbbells, hard triatomics, mixture of hard dumbells of different site-site distances, and n-butane taken as a mixture of conformers. In all the cases prediction within error bars is obtained. (C) 2003 American Institute of Physics.