A general nonorthogonal coupled-electron pair approach based on the intermediate optimization of virtual orbitals is presented. The resulting procedure, similar to the independent electron pair approximation scheme, is developed in the framework of the valence bond (VB) theory, where the effect of the overlap is directly taken into account. Nonorthogonal virtual orbitals optimal for intermolecular correlation effects were determined starting from the self-consistent field for molecular interaction wave function. These were used in the context of a general ab initio variational multistructure VB wave function consisting of double excitations arising from simultaneous single excitations localized on each monomer. The basis set superposition error is excluded in an a priori fashion and geometry relaxation effects are naturally taken into account. As an application example, the equilibrium structure and binding energy of the water dimer system were determined. The equilibrium oxygen-oxygen distance results to be 2.954 Angstrom, in good agreement with the experimental values (2.946 or 2.952 Angstrom) corrected for anharmonicity of the dimer vibrations. The estimated equilibrium interaction energy is -5.02 kcal/mol, thus comparing favorably with the experimental value of -5.44 +/- 0.7 kcal/mol. Taking zero-point vibrational effects into account, the calculated binding enthalpy is -3.22 kcal/mol, in accordance with the experimental estimate of -3.59 +/- 0.5 kcal/mol, determined from measures of thermal conductivity of the vapor. The importance of employing basis sets that include diffuse polarization functions in correlated calculations on hydrogen-bonded systems is confirmed.