In this work, we present a fully automated method for the construction of chemically meaningful sets of non-redundant internal coordinates (also commonly denoted as Z-matrices) from the cartesian coordinates of a molecular system. Particular focus is placed on avoiding ill-definitions of angles and dihedrals due to linear arrangements of atoms, to consistently guarantee a well-defined transformation to cartesian coordinates, even after structural changes. The representations thus obtained are particularly well suited for pathway construction in double-ended methods for transition state search and optimisations with non-linear constraints. Analytical gradients for the transformation between the coordinate systems were derived for the first time, which allows analytical geometry optimizations purely in Z-matrix coordinates. The geometry optimisation was coupled with a Symbolic Algebra package to support arbitrary non-linear constraints in Z-matrix coordinates, while retaining analytical energy gradient conversion. Sample applications are provided for a number of common chemical reactions and illustrative examples where these new algorithms can be used to automatically produce chemically reasonable structure interpolations, or to perform non-linearly constrained optimisations of molecules.