What minimization are you talking about?

In principle, every quantum chemical method is capable of calculating electron densities, or to be precise, probability distribution functions.

The difference between the methods is just the variation of the accuracy.

We just want to calculate the density function of a particular confirmation state (i.e. such as bioactive conformer/co-crystallized ligand). How it can be done using Gaussian software without changing the conformation of co-crystal ligand?

If you want to freeze one ligand moving only the other, you can fixate their coordinates:

https://gaussian.com/opt/ [third tab: options]

My recommendation would be to preoptimize the structure with all other coordinates frozen and then allow the intermediate structure to relax because otherwise you will get structures that may be too far from reality. It would be logical to assume that if one ligand changes its conformer, the other ones will at least elongate or shorten some bonds.

Thanks for your suggestion. I will go through tge website as you have suggested.