The calculation of molecular Hessians for large-scale multiconfiguration self-consistent-field (MCSCF) functions is described. The formalism is based on exponential parametrization of the wave function and symmetric orthonormalization of the molecular orbitals. Extensive use is made of one-index transformations of the molecular integrals, both to construct the gradient vectors that appear in the linear MCSCF response equations, and to perform the multiplication of the trial vectors on the electronic Hessian in the iterative, direct solution of the response equations. No element of the electronic Hessian is ever calculated explicitly, allowing for use of large configuration expansions. Efficient methods are developed for obtaining the solution vectors of the linear response equations. The accuracy of the molecular Hessian is analyzed in terms of the accuracy of these solution vectors. To allow for large basis sets Fock matrices are used to minimize transformations and integrals are recalculated to minimize storage requirements. Integral derivatives are calculated following the McMurchie-Davidson scheme. A simplified algorithm for calculation of derivatives of integrals involving one-center overlap distributions is described. Sample calculations involving several thousand configurations are reported.