I am confused about the relation between the natural and SI units in the special relativity. For example, I have found that I am not able to solve a classical problem of finding the center-of-mass energy of a head on collision of two protons in the SI units. I defined the four vector of proton 1 and 2 to be: $\vec{P_1} = (E_p/c, \vec{p})$ and $\vec{P_2} = (E_p/c, \vec{-p})$
Then I defined the square of the center-of-mass energy to be:
$s = (\vec{P_1} + \vec{P_2})^2 = |\vec{P_1}|^2 + |\vec{P_2}|^2 + 2 \vec{P_1} \vec{P_2} = (E_p^2/c^2 - |\vec{p}|^2) + (E_p^2/c^2 - |\vec{p}|^2) + 2 (E_p^2/c^2 - \vec{p}\vec{p} \cos(90 ^{\circ}))$
$s = 4E_p^2/c^2$
$\sqrt{s} = 2E_p/c$ in the SI units ($\sqrt{s} = 2 E_p$ in the natural units)
However, the equation derived in the SI units seems to be wrong. For example, assuming $E_p$ to be 1 GeV = $10^9 \times 1.602 \times 10^{-19} = 1.602 \times 10^{-10} J$. Then $\sqrt{s} = 2E_p/c = \frac{2 \times 1.602 \times 10^{-10}}{3 \times 10^8} = 1.068 \times 10^{-18}$ J $= 6.7 eV$ which is obviously wrong.
Thank you for your help.