Construction of a valid interaction potential function between quarks is a crucial issue in hadronic physics and also one of the frontier issues. Non-relativistic Breit potential is a common model to describe the interaction between quarks. It is used to successfully calculate the bound states of quarks and quark scatterings. These spur people to improve it. As is well known, the full Breit potential function, which includes the color-Coulomb term, the mass term, the orbit-orbit interaction term, the spin-spin interaction term, the spin-orbit interaction term, the tensor force term, and the constant term, contains singularity factors. How to eliminate the singularity factors is the most urgent task for developing Breit potential model. In this paper, we carry out a replacement method to eliminate the singularity factors in the full Breit quark potential function in coordinate space. Except for the color-Coulomb term and the constant term, remaining terms in the Breit quark potential function are all reconstructed. The replacement of (r) 3 e-r/8 is applied to the mass term and the spin-spin interaction term. The replacement of 1/r (1-e-r)/r is applied to the obit-obit interaction term. The replacement of 1/r3[1-(1+r)e-r]/r3 is applied to the spin-obit interaction term and the tensor force term. We calculate mass splits of heavy mesons and quarkonium species by using the reconstructed potential function and test the validity of the reconstructed potential function. The screening mass used in the calculations is not a simple constant but a variable relating to the quark mass mi and mj. It is found that the simple screening-mass expression cannot give the accurate value of B-meson mass, although it may give the mass splits of light mesons. However, the calculated results of the mass splits of the light mesons -, the heavy mesons, c-J/, b-(1s), c0-c2, etc., are highly consistent with the experimental data only when the screening mass is taken to be the Laurent series, =c-3(a+0.512)-3+ c-2(a+0.512)2 +c-1(a+0.512)-1+c0+c1(a+0.512) with respect to the average quark mass a=(mi+mj)/2. In this case, the mass accuracy of other mesons, especially the six D mesons, is improved significantly. Our calculated results indicate that a valid quark potential model, which gives not only the mass values of light mesons accurately but also the mass splits of heavy quarkonium species, is thus constructed in this paper.