Modeling in biological systems is important at every level from the molecular, to the cellular, to the tissue level. In this paper we discuss the following problem in molecular modeling: given a three-dimensional conformation of a molecule, how do we automatially compute the conformations of the molecule that satisfy certain spatial constraints, that is, certain "feature" atoms of the molecule are in user-specified positions? This task is important in the analysis of receptor-ligand interactions and in other applications such as drug design and protein folding. Using an analogy between robots and molecules, we use the term inverse kinematics to describe the above conformational problems. To solve these problems, we first derive a system of polynomial equations. Then, we adopt a technique based on the Groebner basis from algebraic geometry and develop a novel subdivision algorithm to approximate the real solutions. The approximated solutions can then be used as the starting conformations for existing (heuristic) energy minimization procedures that try to satisfy the target positions of feature atoms and reduce the overall energy of the conformation. To our knowledge, this is the first time that a rigorous algebraic methodology has been used to approximate molecular inverse kinematics solutions.