A new approach to calculate nuclear magnetic resonance chemical shifts has been implemented at the semiempirical modified neglect of diatomic overlap level using gauge-including atomic orbitals. The perturbed density matrix with respect to the magnetic field is obtained by the diagonalization of the complex Fock matrix using the divide and conquer (DC) method, instead of by solving the computationally expensive coupled perturbed Hartree-Fock equations. Adopting the Patchkovskii and Thiel parameters [S. Patchkovskii and W. Thiel J. Comput. Chem. 20, 1220 (1999)], we were able to reproduce their results for small organic molecules. The errors introduced by DC method are negligible, as shown by the calculations on a series of polyalaine structures. Test calculations on proteins have demonstrated that our approach makes it possible to calculate chemical shifts routinely on systems with hundreds of atoms with good accuracy. (C) 2004 American Institute of Physics.