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

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.
• No worries about that, I'm just telling you to help. Before we do a deep dive into the level, I suggest reading DFT Functional Selection Criteria, where I and others have gathered quite a bit on the matter. Tl;dr: It's complicated and depends on what you are looking for. What you most definitely will need is some dispersion correction. But I guess we'll burn that bridge once we've crossed it. In its core, the question is good, so let's hope someone can provide an answer. – Martin - マーチン Feb 12 '20 at 17:39
• Also (loosely) related: How to calculate Lennard-Jones potential with quantum mechanical methods, How do water models fit into force field methods? I was pretty sure there are more related questions, but I just can't find them. I recommend browsing around with the search. – Martin - マーチン Feb 12 '20 at 17:57
• Thanks. I saw this reference, but my question, at first, would be how to build a coherent simulation in Gaussian to obtain the energy curve! "what should i take into account in this kind of simulation?" I mean. As you said, two things are important: level of theory and dispersion correction. – Emerson P L Feb 12 '20 at 18:03
• Is it for fun or for production? You are aware that real dispersion forces are multi-body interactions and DFT does not calculate them, right? – Greg Feb 16 '20 at 16:52

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:

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.

• I've probably missed some of the various "fit force field from QM" papers since it's a bit outside my expertise. Suggestions welcome. – Geoff Hutchison Feb 13 '20 at 16:16
• Even if you do your own fitting, I suggest at least reading these papers, since they've all spent time considering requirements. – Geoff Hutchison Feb 13 '20 at 16:17
• Thanks a lot Geoff Hutchison! I will certainly read the articles you recommended. Your comments are an excellent starting point. – Emerson P L Feb 13 '20 at 19:44
• @EmersonPL - Because HF doesn't fully include electron correlation, it does not properly handle dispersion either. There are empirical dispersion corrections for HF as well, e.g. HF-3c, HFD, .. depending on the program you use. – Geoff Hutchison Feb 23 '20 at 0:45
• If you have the time / computational resources, I'd recommend an RI-MP2 method, eg with cc-VDZ or def2-SVP since it would also include better treatment of dispersion. – Geoff Hutchison Feb 23 '20 at 0:46