The coupled cluster method with single and double excitations has been formulated in a basis set independent language of first quantization. In this formulation the excitation operators are represented in terms of one- and two-electron cluster functions satisfying a set of integrodifferential equations and the strong orthogonality conditions. These equations are solved iteratively by minimizing appropriate Hylleraas-type functionals. During the iteration process correlation energies of up to fourth order in the Moller-Plesset perturbation operator are extracted. A slight modification of the coupled cluster equations leads to an explicitly correlated formulation of the configuration interaction theory. The method was tested in applications to two- and four-electron systems: He, Li+, H-2, Be, Li-, and LiH. The two- electron cluster functions were expanded using explicitly correlated Gaussian geminal bases optimized in the lowest order of perturbation theory. Most of the correlation energies computed at various levels of the coupled cluster and perturbation theory represent the most accurate values to date.