Dynamics of proteins are crucial for its ability to perform its functions, but are often overlooked in the plethora of protein function prediction methods. However, the need for faster screening over sequence-structure-function spaces is ever growing due to the increase in number of known peptide sequences and structures without function annotations. Machine learning methods are commonly implemented to meet the latter need to moderate success, with graph neural networks (GNN) gaining traction in recent years due to its rather intuitive mapping to protein structures. On the other hand, encoding the dynamics of proteins is yet a critical outstanding issue. In our attempt to close this gap, we encode the dynamical information of modal shapes and frequencies, along with the protein contact map, into multigraphs to aid graph neural networks in their prediction. This reduction of mode shapes into graphs is done by calculating the correlation between pairwise residue movements, then weight-averaging the correlation across multiple mode shapes, taking into consideration the equipartition theorem. The anisotropic network model (ANM) and the torsional network model (TNM) are implemented for modal analysis, and the two network models are compared within our framework. The method is tested on various datasets with different sequence similarity cutoffs, and the results compared with existing methods. Furthermore, we are able to extract from our model the dynamically-critical residues, such as hydrogen-binding and ligand-binding sites, for each of the protein's functions. As such, the proposed model provides insight to critical mutation points, facilitating further investigation and understanding of the protein.