Hydrogen bonding of DNA bases was investigated by reliable nonempirical ab initio calculations. Gradient optimization was carried out on 30 DNA base pairs using the Hartree-Fock (HF) approximation and the 6-31G** basis set of atomic orbitals. The optimizations were performed within C-s symmetry. However, the harmonic vibrational analysis indicates that 13 of the studied base pairs are intrinsically nonplanar. Interaction energies of base pairs were then evaluated at the planar optimized geometries with inclusion of the electron correlation energy using the second-order Moller-Plesset (MP2) method. The stabilization energies of the studied base pairs range from -24 to -9 kcal/mol, and the calculated gas phase interaction enthalpies agree well (within 2 kcal/mol) with the available experimental values. The binding energies and molecular structures of the base pairs are not determined solely by the hydrogen bonds, but they are also strongly influenced by the polarity of the monomers and by a wide variety of secondary long-range electrostatic interactions that also involve the hydrogen atoms bonded to ring carbon atoms. The stabilization of the base pairs is dominated by the Hartree-Fock interaction energy. This result confirms that the stability of the base pairs originates in the electrostatic interactions. For weakly bonded base pairs, the correlation interaction energy amounts to as much as 30-40% of the stabilization. For some other bask pairs, however, repulsive correlation interaction energy was found. The latter fact is explained as a result of a reduction of the electrostatic attraction upon inclusion of the electron correlation. The empirical London dispersion energy does not reproduce the correlation interaction energy. For the sake of comparison, results of a first gradient optimization for a DNA base pair at a correlated level (CC base pair, MP2/6-31G** level) are reported. In addition, the ability of the economical density functional theory (DFT) method to reproduce the ab initio data was investigated. The DFT method with present functionals is not suitable to consistently study the whole range of the DNA base interactions. However, it gives good estimates of interaction energies at the reference HF/6-31G** geometries.