Cannabinoids are a group of chemical compounds that have been used for thousands of years due to their psychoactive function and systemic physiological effects. There are at least two types of cannabinoid receptors, CB1 and CB2, which belong to the G protein-coupled receptor superfamily and can trigger different signaling pathways to exert their physiological functions. In this study, several representative agonists and antagonists of both CB1 and CB2 were systematically studied to predict their binding affinities and selectivity against both cannabinoid receptors using a set of hierarchical molecular modeling and simulation techniques, including homology modeling, molecular docking, molecular dynamics (MD) simulations and end point binding free energy calculations using the molecular mechanics/Poisson-Boltzmann surface area-WSAS (MM-PBSA-WSAS) method, and molecular mechanics/generalized Born surface area (MM-GBSA) free energy decomposition. Encouragingly, the calculated binding free energies correlated very well with the experimental values and the correlation coefficient square (), 0.60, was much higher than that of an efficient but less accurate docking scoring function ( = 0.37). The hotspot residues for CB1 and CB2 in both active and inactive conformations were identified MM-GBSA free energy decomposition analysis. The comparisons of binding free energies, ligand-receptor interaction patterns, and hotspot residues among the four systems, namely, agonist-bound CB1, agonist-bound CB2, antagonist-bound CB1, and antagonist-bound CB2, enabled us to investigate and identify distinct binding features of these four systems, with which one can rationally design potent, selective, and function-specific modulators for the cannabinoid receptors.