We present detailed methodology and results for the approximate calculation of the anharmonic vibrational wave functions for a given structure of a protein. The ground state and the fundamental vibrationally excited states of hydrated BPTI are computed with an approach that is of very good accuracy for low-lying states. The eigenfunctions are used to predict quantum properties such as vibrational excitation frequencies and IR intensities as well as atomic mean square displacements at T = 0 K. The method treats diagonal anharmonic effects exactly up to fourth order in normal mode coordinates, while using a mean-field approximation, the vibrational self-consistent field (SCF) approximation for the mode-mode couplings. The inclusion of the diagonal effects (exact in fourth order) makes the potentials stiffer than quadratic. When the mode-mode coupling is included as a mean field, the system becomes even stiffer. These effects produce significant modifications of such observables as the vibrational spectrum and Debye-Waller factors. The results presented here should be considered essentially exact except for potentially important approximations : (a) that the potential energy used bears some resemblance to reality, (b) that the system is restricted to a single minimum in the potential energy surface, and (c) that the quartic expansion in normal modes, when truncated to fourth order, provides a good description of the potential. These issues will be addressed in the main text. The results show that the vibrational absorption spectrum is strongly affected by anharmonic and mode-mode coupling effects. Deviations from the corresponding harmonic absorption intensities are very large. Unlike our previous study that found only weak mode-mode coupling effects for the lowest 100 modes, the SCF corrections calculated here for the intermediate frequency modes are seen to be very significant. The results have important implications for the vibrational spectroscopy and other properties of proteins at low temperatures.