Computation of the Fock matrix is currently the limiting factor in the application of Hartree-Fock and hybrid Hartree-Fock/density functional theories to larger systems. Computation of the Fock matrix is dominated by calculation of the Coulomb and exchange matrices. With conventional Gaussian-based methods, computation of the Fock matrix typically scales as similar to N-2.7, where N is the number of basis functions. A hierarchical multipole method is developed for fast computation of the Coulomb matrix. This method, together with a recently described approach to computing the Hartree-Fock exchange matrix of insulators [J. Chem. Phys. 105, 2726 (1900)], leads to a linear scaling algorithm for calculation of the Fock matrix. Linear scaling computation the Fock matrix is demonstrated for a sequence of water clusters at the restricted Hartree-Fock/3-21G level of theory, and corresponding accuracies in converged total energies are shown to be comparable with those obtained from standard quantum chemistry programs. Restricted Hartree-Fock/3-21G calculations on several proteins of current interest are documented, including endothelin, charybdotoxin, and the tetramerization monomer of P53. The P53 calculation, involving 698 atoms and 3836 basis functions, may be the largest Hartree-Fock calculation to date. The electrostatic potentials of charybdotoxin and the tetramerization monomer of P53 are visualized and the results are related to molecular function.