Normal coordinate analysis (NCA) requires the diagonalization of a mass-weighted potential energy second derivative matrix (Hessian matrix).  In chemical systems with thousands of atoms, the Hessian matrix is extremely sensitive to small changes in structure, leading to negative eigenvalues (imaginary frequencies).  This problem can be eliminated by performing a molecular dynamics simulation and constructing a trajectory averaged Hessian matrix.  With this method and using sparse matrix methods, normal coordinates and frequqncies have been calculated for polymer particles and crystals with up to 24,000 atoms.  Computation of heat capacities and other thermodynamic quantities that can be computed from g(omega) are described in detail.  Also described are approximate methods for obtaining thermodynamic quantities that do not require full diagonalization of the Hessian matrix.