Question: I would like to construct a model analytical potential function for collision between CH3Cl and Ar. This would be the sum of CH3Cl intramolecular + CH3Cl---Ar intermolecular potential. First I would like to generate the ab Initio data points to be fitted using double many-body expansion or similar functions. I know how to set up PES scans (bond, angle, dihedral etc) using gaussian, Orca and Psi4 but I don't know how I would approach this problem. I would greatly appreciate if anyone would give a detailed procedure on how I can go about solving this and similar problems. Thanks.

Answer: 
We attempted to solve a similar problem when studying the (also highly symmetrical) $\ce{CH4}$ and $\ce{CF4}$ homo- and hetero-dimers. I found it easiest to use internal coordinates and fixed molecular geometries to generate the hypersurface, but we also had large computational resources at our disposal and we might have benefitted from using even more, especially to generate more points at larger intermolecular distance.[Chattoraj, Risthaus, Rubner, Heuer, Grimme, *J. Chem. Phys.*, **142**, 164508 (2015)](http://dx.doi.org/10.1063/1.4919079)

In your particular case, you should be able to use

* the $\ce{Ar}$-$\ce{C}$ distance
* the $\ce{Ar}$-$\ce{C}$-$\ce{Cl}$ (0 to 180°)
* a $\ce{Ar}$-$\ce{C}$-$\ce{Cl}$-$\ce{H}$ angle, wherein you use the $\ce{H}$ closest to the $\ce{Ar}$ (-60 to 60°)

The tricky part remains: picking the right values with which to cover the intervals, which is a bit of trial and error. Some researchers have used "self-healing" approaches, which mean that the configurations with the worst differences between reference and your model are used to refine the model.

If you plan to not use a fixed molecular geometry, my approach is of course not valid.
On the fitting, I would recommend Dr. Nike Dattani's answer. You can also try some other methods such as those using neural networks.

I read your discussion and find that may be you have also some issues with the scanning using ab initio methods.

The intermolecular interaction is often read
$$\Delta V\_{\rm int}=\Delta V\_{\rm int}(R, \theta, \phi),$$
where $R$, $\theta$, and $\phi$ define the position of Ar as if the COM of CH3Cl locates at the origin of the polar coordinate system.
To scan this, the most common way to do the job is to write the Z-matrix,
which defines the relative position of the atoms by given the distance between two, angle between three, and dihedral between four atoms.
In our group, this job is often done by hand, because there are no to much atoms here.
If you need, you can check the manual of your preferred software, like in Molpro, it is [here](https://www.molpro.net/manual/doku.php?id=molecular_geometry#z-matrix_input).

If you can define the geometry of one geometry, you can scan the global PES by changing the parameters in the Z-matrix.

I believe the trick here is to use `dummy centre` to represent the COM. It means the centre is just an indicator for a reference position instead of a real atom. In the link to the Molpro manual, it is the `Q` or `X` atom.

Also notice that if you are scanning the intermolecular PES, you need to minimize the Basis set suerposition error (BSSE) by using [counterpoise correction](https://www.molpro.net/manual/doku.php?id=molecular_geometry#dummy_centres),
where you need calculate the energy using the same basis (remember to include the basis for the dummy molecule) three times for every geometry and get the interaction potential as
$$\Delta V\_{\rm int}=E({\rm CH\_3Cl+Ar})-E({\rm CH\_3Cl})-E({\rm Ar}),$$
where the $E$'s are the energy get from the ab initio package you use.

Note that you will find the terminology are greatly different among packages.
In my limited experience, [CFOUR](http://slater.chemie.uni-mainz.de/cfour/index.php?n=Main.MolecularGeometryInput) will call the first dummy centre (the `Q` or `X` ones) as 'dummy atom' while call the latter one 'ghost atoms'. The difference between the two lies in whether the basis functions of the atom are kept in the ab initio computation.