A method is given for the exact calculation, including all anharmonic effects, of the matrix elements of the true Hamiltonian H of a crystalline solid between the eigenfunctions of any arbitrary harmonic Hamiltonian appropriate to the crystal structure. By the introduction of a variational scaling parameter a, a set of harmonic eigenfunctions can be generated from the eigenfunctions of a single such harmonic Hamiltonian. The ground-state energy W0 of the crystal and the optimum value a0 of a are then determined by a variational calculation which minimizes the expectation value of H between the ground-state harmonic eigenfunctions. The other diagonal matrix elements of H, calculated using the states determined by a0, then give first-order values for the energy of one-phonon and multiphonon states. In addition, the off-diagonal matrix elements can be obtained and used in a perturbation calculation to improve the energies. It is shown that the value of W0 will always be lower than energy obtained from a closely related variational calculation using a wave function constructed out of a product of Gaussian orbitals centered on the lattice sites. The decrease in energy is about 12% of the original kinetic energy, and is divided approximately equally between the kinetic- and potential-energy contributions. Numerical results are given for W0 and for a few phonon energies of a model bcc crystal in which each atom interacts with its first two shells of nearest neighbors through a realistic atomic potential. Possible applications of this theory to rare-gas solids are discussed. © 1966 The American Physical Society.