Optimizations are all about finding the minimum in something. Typically in geometry optimization, it is about finding the minimum in energy. At a minimum the derivative of energy with respect to changing position should be zero, or better put, the jacobian should be positive definite.
Your question is about bulk systems, however, single molecules are problematic enough.
A single molecule by itself can take on many different conformations, each with a different energy. Given a starting configuration, a geometry optimization changes atomic coordinates to minimize energy. In practice, this usually means you find the nearest local minimum.
It is therefore important, and often ignored, to do a conformer search for a single molecule and find the lowest energy conformer prior to doing a geometry optimization. I can't stress the importance of conformer searches enough.
You are interested in a system of many molecules. For a single molecule the actual x,y,z position does not really matter. For a system of them however it does. Each molecule will relax to some conformer, possibly not the best one, and they will relax to certain orientations and centers of mass x,y,z most likely not the best ones, but representing the nearest local minimum in energy from the initial starting guess.
Finding a best geometry for a single molecule is tough, I would say, currently, it is impossible for a bulk system, you find the best one that you can, and live with it.
You can of course generate many initial guesses and take the lowest energy final geometry. You can try simulated annealing techniques, you can try all sorts of numerical methods really, but, finding the global minumum is an unsolved problem for an N-body problem such as a bulk phase of molecules.
Courtesy of Andrew Rosen, this paper has a nice example of zeolite structures depending on initial configuration.