Even when using parametrized semiempirical methods, quantum chemical calculations on molecules containing more than a few hundred atoms become prohibitively expensive due to O (N3) time and memory costs where N is the number of atoms. I implemented methods to allow the CPU time cost of semiempirical methods to scale linearly with system size enabling semiempirical calculations on large biological systems such as proteins and nucleic acids. The cost of forming the initial guess density matrix was reduced by replacing the O (N3) diagonalization of the Huckel Hamiltonian with an approach which uses localized molecular orbitals based on the Lewis dot structure to build the density matrix. The Fock matrix build was reduced from O (N2) to linear scaling in CPU time using atom-atom distance cutoffs. The diagonalization step was replaced by several linear scaling methods described in the literature: conjugate gradient density matrix search (CGDMS), purification of the density matrix (PDM), pseudodiagonalization (PD), and the Chebyshev expansion method (CEM). While in my semiempirical implementation all of these methods demonstrated linear scaling, CGDMS, PDM and PD required about the same amount of CPU time for calculations on water clusters and polyglycine chains but CEM was found to be about three times as expensive as the other methods. However, CGDMS stands out among the other methods by having the added property of enhancing self-consistent field (SCF) convergence in cases where diagonalization has convergence difficulties. Finally, to demonstrate the effectiveness of the linear scaling semiempirical method on a realistic system we performed the first-ever semiempirical geometry optimization using PM3 implemented with CGDMS on a 1226 atom kringle 1 of plasminogen.