Deformation in an aggregate of single crystalline grains, in conventional crystal plasticity finite element methods, is modeled through the collective behavior of individual ones each consisting of multiple elements. Such a procedure often entails tedious preparedness such as grain-level discretization and post-processing to obtain meaningful collective behavior. Virtual Element Method (VEM) which is designed to take care of a multifaceted three-dimensional body as one element, offers a promising alternative for crystal plasticity modeling. In this paper, we implement a phenomenological, rate-dependent crystal plasticity model within the framework of the energy-based VEM. By integrating nested iterative algorithms with automatic differentiation techniques, the numerical procedure efficiently simulates both single-crystal and polycrystalline solids under various loading conditions. With three numerical examples, we demonstrate the ability of VEM based crystal plasticity modeling to predict stress-strain responses, texture evolutions, and deformation mechanisms in face-centered cubic crystals. This numerical procedure shows a promising method of modeling polycrystalline solids composed of a vast amount of grains, and paves a new route of two-scale modeling which bridges microstructures with macroscopic mechanical behavior of solids.