The method of the bond function basis set combined with the counterpoise procedure is studied in detail by the complete fourth-order Moller-Plesset perturbation (MP4) theory, following from a recent communication report [J. Chem. Phys. 98, 2481 (1993)]. This method is applied to calculate molecular dissociation energies D-e as well as equilibrium bond distances r(e) and harmonic frequencies omega(e) of a number of diatomic molecules (N-2, O-2, F-2, Cl-2, HF, HCl, and CO) and the results are compared with those from other methods, without either counterpoise procedure or bond functions or both. The usefulness of the method is shown by the results for all the molecules using a moderately polarized basis set (2p1d for H atom and 2d1f for heavy atoms) augmented with the universal bond functions 3s3p2d. The method has consistently recovered 98%-99% of the experimental values for D-e, compared to as low as 40% without bond functions. The effect of bond functions is less significant on the predictions of r(e) and omega(e), due primarily to the inadequacy of the MP4 theory, but our method is still shown to be favored over the other methods. The electric dipole moments of the polar molecules (HF, HCl, and CO) are also examined and it is found that the use of bond functions results in a significant improvement of the dipole values. Detailed discussions are given to explain the need for bond functions and the counterpoise procedure. The high linear independence with nucleus-centered basis functions is explained to be responsible for the efficiency of bond functions. The counterpoise procedure is logically justified from the conventional noncounterpoise procedure. Potential problems and limitations associated with the proposed method are also discussed.