The three-photon decay rate of J=ψ is studied using two N f ¼ 2 twisted-mass gauge ensembles with lattice spacings a ≃ 0.085 and 0.067 fm. Using a new method, only the correlation functions directly related to the physical decay width are computed with all polarizations of the initial and final states summed over. The final result can be obtained after a naive extrapolation to the continuum limit. To be specific, the result for such rare decay is given as BðJ=ψ → 3γÞ ¼ ð2.13 AE 0.14 AE 0.29Þ × 10 −5 , where the first error is statistical and the second is an estimate of the systematics. We also propose a method to analyze the Dalitz plot of the corresponding process based on the lattice data which can provide direct information for the relevant experiments.