How to obtain curve energy and non-bonded interaction parameters (e.g. Lennard Jones parameters) with Gaussian?

Question

I am parameterizing a force field and would need to obtain non-bonded interaction parameters of atoms (e.g., Lennard-Jones parameters).

My primary idea is to obtain a potential energy curve using the Gaussian 09 software and fit the data on that curve to an equation that has non-bonded interaction parameters of atoms (e.g., Lennard-Jones curve).

What I thought was to put two identical molecules and vary the distance between two atoms of these molecules (e.g., using the scan function on redundant coordinator editor) and, therefore, get the energy curve. What I was trying to do was a simple example to get interaction parameters for the oxygen atoms of two water molecules.

I was taking a test with Gaussian at the theory level B3LYP/6-311G++(d,p). Here is the input file:

%nprocshared=4 
# opt=modredundant b3lyp/6-311++g(d,p) geom=connectivity

Title Card Required

0 1
 O                 -3.09468841   -0.10392610    0.00000000
 H                 -2.13468841   -0.10392610    0.00000000
 H                 -3.41514300    0.80100974    0.00000000
 O                  3.39491938    0.10392610    0.00000000
 H                  4.35491938    0.10392610    0.00000000
 H                  3.07446479    1.00886193    0.00000000

 1 2 1.0 3 1.0
 2
 3
 4 5 1.0 6 1.0
 5
 6

 B 1 4 S 10 0.200000

Something I noticed is that the level of theory used (obviously) influences the calculations results. In addition, the initial distance of the water molecules also influences whether the simulation will fail or not (in this example, it works).

Could you kindly tell me if this approach makes sense? Or should I use another software, another type of calculation? Do you have any suggestions on how I could get an energy curve by varying the distance of two non-bonded atoms in the Gaussian? Examples would be appreciated.

With Gaussian and a similar approach I was able to obtain parameters such as bond stretching, angle deformation and dihedral torsion. However, it would be lacking to obtain the non-bonded interaction parameters of atoms.

Things I should consider:

1. Do a rigid scan instead a relaxed scan (like the code presented).
2. I need a dispersion-corrected density functional methods.

Answer 1

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.