A group theoretical method is described for calculating the normal modes of icosahedral systems such as viruses. The use of symmetry reduces the size of the matrices that have to be diagonalized from 60Nx60N to 5Nx5N (where N is the number of degrees of freedom of every subunit) while preserving complete accuracy. The method includes algorithms to calculate the normal modes, the atomic fluctuations and cross-correlations, and the projections of normal modes on vectors of interest. The correctness and accuracy of the method are verified using a model system, consisting of 60 Dialanine peptides arranged in a nonbonded icosahedral complex. The effects of using reduced basis sets in the normal mode calculations are examined. Reduced basis sets, especially those consisting of dihedral and/or bond angles, are shown to have relatively small effects on the frequencies, spatial fluctuations, and directions of the normal mode displacements. The current implementation allows accurate reduced basis normal mode calculations on icosahedral virus molecules with moderately powerful computers.