A strategy is proposed to determine atomic wave functions by the configuration-interaction method (CI). A large nonoptimized set of Slater-type orbitals (STO's) is contracted into an orthogonal set of smaller dimension whose elements are natural orbitals (NO's) of a CI wave function which includes the Hartree-Fock (HF) configuration and, out of all possible double excitations, only those such that the HF orbitals are also NO's. Conceptual and practical arguments which justify such a procedure are given. A nonrelativistic calculation of the ground state of Be is made with an STO basis used earlier by Watson, but including s, p, and d orbitals only. The final 180-term CI expansion (1492 Slater determinants) gives an energy E=-14.664193 atomic units (a.u.) (Be), which is about 0.00221 a.u. higher than a nonrelativistic estimate for this state, and the pair energies ∊(1s, 1s)=-0.040869 a.u., ∊(2s, 2s)=-0.045104 a.u., and ∊(1s, 2s)=-0.005240 a.u., are in excellent agreement with those found by Kelly by means of a many-body perturbation calculation, and also by Byron and Joachain, who used a variation-perturabtion method. Quantitative comparisons with all previous work on Be are made. Single excitations occur with large eigenvector components, but they become vanishingly small when NO's of the final wave function are used as the basis. The eigenvector components of quadruply excited configurations are closely related to those arising from a separated-electron-pair wave function without the strong orthogonality condition. The spd energy limit is estimated to be Espd=-14.66453 a.u. From considerations not involving relativistic corrections, the "exact" nonrelativistic energy is estimated to be lower than E=-14.66639 a.u. by no more than 0.0003 a.u. The "exact" pair energies are estimated to be ∊(1s, 1s)=-0.04261 a.u., ∊(2s, 2s)=-0.04550 a.u., and ∊(1s, 2s)=-0.00530 a.u. Details on all aspects of the calculation are given. Further work on several states of the first-row atoms is in progress.