Your software library uses a gradient descent method for inverse kinematics. The way this works is that it moves the Cartesian position of the end-effector towards the goal in small steps then gets the joint angles for each step with the relation:
$$
\Delta \theta= J^{-1}\Delta x
$$
This method has a number of nice properties, but as you found you only get one solution.
An analytical inverse kinematics solution is what you want. In this method, you algebraically determine all possible solutions. The only catch is if you are under-constrained, (like in your case), where you have more joints than degrees of freedom. (You have 6 joints, but only want a 3 DOF position). In this case, you pick some joints to "fix" and do the math as if that joint is rigid. But in reality, you typically discretize that joint's range, and get a set of solutions for every fixed position. (You can see how the curse of dimenionality will get you here.)
I would recommend investigating ikfast from OpenRave, which is an analytical IK solver. This has a number of different IK types for different DOFs. "Translation3D" is what you want.
But the way to trick your gradient decent IK solver into giving you different solutions is to give it different initial conditions. You should start from drastically different arm configurations. (Like joints flipped 180 degrees from each other).
As for the not needing rotation: The gradient descent IK method can ignore rotation quite easily. You would simply only use the top 3 rows of the Jacobian, and give it only the end-effector position (not rotation too).
