A method is developed to analytically evaluate energy derivatives for extended systems. Linear dependence among basis functions, which almost always occurs in extended systems and brings instability to the coupled-perturbed equations, is automatically eliminated in this method. The remaining independent basis functions are transformed into semiorthogonal orbitals. The derivatives of the orbitals and the overlap matrix over them are obtained via a set of coupled-perturbed equations, similar to those of the coupled-perturbed Hartree-Fock (CPHF) equations which are used to calculate the derivatives of the Hartree-Fock (HF) orbitals and the orbital energies. By introducing symmetrized coordinates. these coupled-perturbed equations can be easily solved. Explicit expressions for calculating gradients and Hessians of the HF energy for extended systems are given. With this method, we can calculate energy derivatives with respect to displacements of the nuclei, including those which break the translational symmetry. Therefore, the method not only provides an efficient and accurate approach to calculate energy derivatives of any order, but also enables us to determine the force constants for individual nuclei, the interatomic force constants, and phonon dispersion curves in the whole Brillouin zone. With this method, the computational cost to calculate phonon spectrum with k not equal 0 in the Brillouin zone is the same as that needed for the spectrum at k = 0.