How to obtain curve energy and non-bonded interaction parameters (e.g. Lennard Jones parameters) with Gaussian?
To start off, you definitely need to use dispersion-corrected density functional methods if you want to obtain reasonable non-bonded interaction parameters. I've seen many cases where B3LYP (with no dispersion) yields a fully-repulsive potential energy scan.

I also agree with your assertion about doing a rigid scan. In principal, a relaxed scan is a good idea, but consider that most force fields separate out the non-bonded interactions. Thus, for parameterizing, you don't really want to change the bond lengths, angles, etc. from a relaxed scan.

There are a few ways to get the parameters. As mentioned above, there's already a similar question with tips on curve fitting.

Multiple papers have also been published on fitting force field parameters from quantum chemical results, e.g: 


QubeKit J. Chem. Inf. Model. 2019, 59, 4, 1366-1381 - code at GitHub
QuickFF J. Comput. Chem. 2015, 36, 1015– 1027 - code at GitHub
ForceBalance J. Chem. Theory Compute. 2013, 9, 1, 452-460 and J. Phys. Chem. Lett. 2014, 5, 11, 1885-1891 - code at GitHub
ForceFit J. Comput. Chem., 2010 31: 2307-2316
Parfit J. Chem. Inf. Model. 2017, 57, 3, 391-396 - code at GitHub


Update there's a paper on this very specific problem:

J. Chem. Theory Comput. 2020, 16, 2, 1115-1127


  Here, we utilize a previously described Minimal Basis Iterative Stockholder (MBIS) method to carry out an atoms-in-molecules partitioning of electron densities. Information from these atomic densities is then mapped to Lennard-Jones parameters using a set of mapping parameters much smaller than the typical number of atom types in a force field. This approach is advantageous in two ways: it eliminates atom types by allowing each atom to have unique Lennard-Jones parameters, and it greatly reduces the number of parameters to be optimized.