We have previously developed a saddle point optimization algorithm that efficiently identifies the leftmost eigenvector of the Hessian matrix using iterative diagonalization. Furthermore, this algorithm uses the information obtained by the iterative diagonalization procedure to construct an approximate Hessian that is exact in the subspace searched by the eigensolver. Our method can be compared to that of Sharada, Bell, and Head-Gordon, with the key difference being that all curvature information obtained by the eigensolver is incorporated into the approximate Hessian, not just the lowest curvature mode. This results in a much more accurate approximate Hessian matrix, and therefore a substantially more efficient optimization algorithm. However, as our method was originally designed only for optimization in a Cartesian coordinate basis, it was not suitable for optimization of molecular saddle point geometries. In this work, we extend our iterative diagonalization and approximate Hessian construction methods to optimization in a basis of redundant internal coordinates, potentially with added constraints. Our new method diagonalizes the Hessian and constructs the approximate Hessian directly in an internal coordinate basis. This is in contrast to the method of Sharada et al., which diagonalizes the Hessian and constructs the approximate Hessian in a Cartesian basis before transforming it into an internal coordinate basis.