Structure of fluids is suitably characterized by distribution functions from which the most important is the pair correlation function. Theoretical approaches to get the pair distribution function are based mainly on the solution of the Ornstein-Zernike (OZ) integral equation. In this paper, the OZ equation for molecular fluids is modified to yield the average correlation function for systems of convex molecules. In our approach we employed the previously proposed method to separate the shape effect of molecular cores from that due to the variable surface-surface distances among three pairs of convex cores. The effect of nonspherical shape of hard cores in the convolution integral is expressed through the derivative with respect to three surface-surface distances of the expression for the hard convex body third virial coefficient. For simple fluids (with the pointwise cores) the derived expression reduces to the standard OZ equation. The modified OZ equation is solved numerically for the Percus-Yevick-type closure and the average correlation functions in the systems of hard spherocylinders with l/sigma =0.4, 0,6 and 1 were determined. The obtained dependencies of the average correlation functions on the reduced distances calculated from the modified OZ equation agree well with the simulation data for the above systems at relatively high densities.