Using Lie algebraic method, we obtain the fitting coefficients of an effective Hamiltonian operator, which conveniently describes vibration spectra of symmetrical linear tetratomic molecules, including both stretching and bending modes. With the classical limit of the Hamiltonian and the kinetic energy equal to zero, we get the potential energy surface of the linear symmetry tetratomic molecule. For an example, we use the method to calculate the potential energy surface of C2D2. We fit 41 data with 7 coefficients; the data are from Iachello et al. [J. Mol. Spectrosc. 149 (1991) 132]. The fitting rms is 6.68cm-1. At last, we use the potential energy surface to get the force constants and the dissociation energy. The method can be applied to a number of linear symmetry tetratomic molecules.