Abstract: The nuclear shell model has been very successful in describing various properties of nuclei, especially in the neighborhood of the closed shells, where the configuration space is usually small but good enough for the construction of the nuclear wave function. However, in the heavy deformed nuclear region, the configuration space is extremely huge. This makes it almost impossible to perform full shell model calculations. To avoid such difficulty, various approximated shell model methods have been developed. Among them, the variation after projection (VAP) method is an important one, in which the nuclear wave functions with good quantum numbers are sufficiently optimized. However, when the VAP method is applied to heavy deformed nuclei, the computational burden is still quite heavy. This is because there are a large number of matrix elements involving the angular momentum projection. All these projected matrix elements should be calculated by integrating the corresponding rotational matrix elements over the three Euler angles. Hence, accelerating the VAP calculation is crucial in extending the VAP application to heavy deformed nuclei, which is the purpose of the present paper. By analyzing the VAP code, it turns out that most of the computational cost is spent on evaluating the projected matrix elements. When the configuration space is huge, such computational cost for the projected matrix elements becomes extremely heavy. Fortunately, all the rotational matrix elements are independent of one another. This makes it very convenient in parallelizing the calculation of these rotational matrix elements. In general, the VAP calculations were accelerated on the traditional CPU platforms. But in recent years, with the rapid developments of the graphics processing unit (GPU), GPU becomes the first choice of modern high-performance computing. Here, by adopting the OpenACC parallel programming directives, the VAP code from the traditional CPU platform to the high-performance GPU platform has been successfully migrated. The calculations of a large amount of the independent rotation matrix elements at each mesh point in the integration of the angular momentum projection were parallelized on a GPU card. It is confirmed that the results calculated by the migrated VAP code are exactly the same as the ones calculated by the original OpenMP parallelized VAP code. Most important of all, the speed of the GPU accelerated VAP code can be several times faster than that of the original OpenMP parallelized VAP code, when they run in the same machine. Using the GPU accelerated VAP code, spectrum of the ground state band of the heavy nucleus ¹⁷⁸Hf have been calculated for the first time. This opens the door of the VAP application to the heavy deformed rare-earth nuclei.