How many coordinates should be used for the potential energy surface of a non-rotating linear molecule? (3N-6) or (3N-5)? Take, for example, $\ce{CO2}$.  It is pretty clear that in the absence of rotation, the energy of this system depends only on the two bond lengths and the bond angle.  

So if I wanted to build a potential energy surface for this molecule that was valid near and far from equilibrium, it should depend on these coordinates.  But how do I rectify this with the fact that $\ce{CO2}$ has four vibrational modes? As we learn in undergraduate chemistry, it has two stretches and two degenerate bends.  

If I build a Hamiltonian in the three valence coordinates, could it possibly describe the ground state energy correctly?  That is, would it contain the zero-point energy for all four modes?  Would it be able to find the right low-lying excited states of the bending mode(s)?

EDIT: 
Wildcat has brought up the point of symmetry and point groups below, and perhaps that can help get to the heart of my question.  When we say that $\ce{CO2}$ is a linear molecule, that means that the lowest energy configuration has a bond angle of $\pi$.  But in no way do I want to restrict the symmetry of the molecule.  I want to work with a potential surface that is valid for symmetric and asymmetric geometries.

But here is where I get confused.  Let's say I did restrict it to be linear.  I would then have a two-dimensional potential, a function of the two bond lengths.  Of the seven normal modes, I only keep two: the symmetric stretch and anti-symmetric stretch.  The ground state energy would be the sum of the zero-point energy of these two modes.   So therefore the number of modes is equal to the number of coordinates I used to write my Hamiltonian.

If I further constrict the symmetry, and impose the $D_{\infty h}$ point group, then the potential is now a function of one parameter.  I will now only have one vibrational mode, the symmetric stretch, and the ground state energy is the zero-point energy of this mode.  So, again, the number of modes is equal to the number of coordinates I use to write my Hamiltonian.

But what if the symmetry were not restricted in any way?  I can write the potential as a function three internal coordinates, be they bond lengths and angle, othogonormal rectilinear coordinates, jacobi coordinates, etc.  This potential, and the Hamiltonian I construct with it, should find the same modes described above, as well as the bending modes.  But will it find both bending modes?  I don't see why it should.  Based on the two reduced-dimensional examples above, I would expect the number of modes to equal the number of coordinates.  I would expect therefore, not to find the right ground-state energy, etc.

Specifically, I am concerned about the accuracy of dynamical simulations performed using this three-dimensional potential.

Likewise, what about a tetratomic liner molecule like $\ce{C2H2}$? I would use three bond lengths, two bond angles, and one dihedral angle, for a total of 6 coordinates - but it has 7 vibrational modes.  Would I model the dynamics properly?
The Cartesian Hessian ($\mathbf{H}_C$) can be obtained by from the Hessian computed in internal coordinates ($\mathbf{H}_i$) through a non-linear transformation
$$\mathbf{H}_C = \mathbf{B}^T\mathbf{H}_i\mathbf{B}$$
where $B_{ij}=\frac{\partial R_i}{x_j}$ and $R_i$ and $x_i$ is an internal and Cartesian coordinate, respectively.

In the case of the degenerate bending mode you can obtain the components of $\mathbf{H}_C$ that you need from $\frac{\partial^2 E}{\partial \theta^2}$, where $\theta$ is the bond angle, through suitable choices of $B_{ij}$'s.  

So the Hessian generated in the three internal coordinates has all the information you need.