Glycogen synthase kinase-3β (GSK3β) is associated with various key biological processes, and it has been considered as a critical therapeutic target for the treatment of many diseases. However, it is a big challenge to develop ATP-competition GSK3β inhibitors because of the high sequence homology with other kinases. In this work, a novel parallel virtual screening strategy based on multiple GSK3β protein structures, integrating molecular docking, complex-based pharmacophore, and naive Bayesian classification, was developed to screen a large chemical database, the 50 compounds with top-scores then underwent a luminescent kinase assay, which led to the discovery of two GSK3β inhibitor hits. The high screening enrichment rate indicates the reliability and practicability of the integrated protocol. Finally, molecular docking and molecular dynamics simulation were employed to investigate the binding modes of the GSK3β inhibitors, and some "hot residues" critical to GSK3β affinity were highlighted. The present study may provide some valuable guidance for the development of novel GSK3β inhibitors.