perfectly within a tiny fraction of a second. Therefore, in the past few decades, many graph theory based biomolecular models, including normal mode analysis (NMA) [@Go:1983; @Tasumi:1982; @Brooks:1983; @Levitt:1985], elastic network model (ENM) [@Bahar:1997; @Bahar:1998; @Atilgan:2001; @Hinsen:1998; @Tama:2001; @LiGH:2002] become very popular for understanding protein flexibility and long time dynamics. In these models, the diagonalization of the interaction matrix or Hamiltonian of a protein is a required procedure to obtain protein eigenmodes and associated eigenvalues. The low order eigenmodes can be interpreted as the slow motions of the protein around the equilibrium state and the Moore-Penrose pseudo-inverse matrix can be used to predict the protein thermal factors, or B-factors. However, NMA interaction potentials are quite complicated. Tirion simplified its complexity by retaining only the harmonic potential for elasticity, which is the dominant term in the MD Hamiltonian [@Tirion:1996]. Network theory [@Flory:1976] has had a considerable impact in protein flexibility analysis. The combination of elasticity and coarse-grained network gives rise to elastic network model (ENM) [@Hinsen:1998]. In this spirit, Gaussian network model (GNM) [@Bahar:1997; @Bahar:1998; @QCui:2010] and anisotropic network model (ANM) [@Atilgan:2001] have been proposed. Yang et al. [@LWYang:2008] have shown that the GNM is about one order more efficient than 