Gaussian scan function help for constructing input file

Question

I have the cartesian coordinates for two molecules that are in a stacked position and I want to be able to calculate the potential energy at 0.2 angstrom steps from when the top molecule is 3.0 angstroms to 7.0 angstroms away in the z direction. My professor told me gaussian has a scan function to do this but I'm having trouble figuring some things out. I have the basis set/ correlation methods I need to use and the cartesian coordinates for all the atoms when they are arranged 3.0 angstroms apart. I tried following this site's advice;

but I'm not sure I constructed the input file correctly. Here is what I have

.

From my understanding, this should just move up all of the atoms of the top molecule by 0.2 angstroms in the Z direction.

Questions:

Do I just throw this stuff anywhere into the input file? What about the original cartesian coordinate array? Surely it has to stay since the "step" array calls back to atoms 101-112. Do I have to change any other commands? The only parameters set are the resource limits and "#P BP86/6-311++G**"

Update:

I think I may have been unclear. I have already optimized my two molecules. One is a nanotube and the other is an arene and my task was to basically generate this plot. I'm trying to find the distance apart where the two molecules are stabilized. Before, I just ran the "#P BP86/6-311++G**"" with my arene @ 3.0 angstrom apart and made new input files changing the arene z coordinates in 0.5 angstrom increments but I can't capture that potential energy well. I decided to try smaller increments to increase the resolution but I have quite a few more situations to test so I wanted to just have Gaussian change the increments and recalculate the potential energy.

Thankfully a full relaxed scan is not necessary in most cases. Oftentimes a nice little rigid scan will give you a decent result. Have a look at the scan keyword. The only downside to this approach is, that you need a z-matrix to make it work. And putting your brain to that task is also a bit tedious.

That's what I initially came across but with ~100 atoms, I just let Chemcraft generate the z matrix from my cartesian coordinates. The thing is, I'm not sure quite how to edit the z-matrix (relevant part below) with the scan commands.

Specifically I want to keep all bond angles the same but move atoms 101-112 in 0.2 A steps in the z direction, and calculate the potential energy.

I guess my question is now if this is a plausible approach. I'll try to look up how to construct the proper z-matrix if that is indeed the best way to do things.

Thanks @Martin

Answer 1

This is (probably) not going to work the way you intend it to work. Constructing a meaningful modredundant calculation is quite a difficult task, as you sometimes cannot anticipate what the program actually does. Using this interface you will perform a relaxed scan, i.e. you will step one variable, fix the values of the scanned coordinates, and run a full optimisation at this point. This allows you to systematically, brute force a potential energy surface. In you case you specify 12 coordinates to scan with 20 points each, resulting in a whopping $N=\frac{20!}{12!}\approx5\times10^8$ optimisations.

Here are a few pointers when constructing a modredundant input:

• don't scan more than two parameters at a time, that roughly equates to a 3D graph of a PES
• don't scan xyz coordinates, scan bond lengths and angles instead; scanning xyz coordinates results in fixing the coordinate system in parts disallowing the algorithm to rotate the molecule/ use coordination transformations to speed up the calculation
• keep constant what you need not optimised; fixing things will result in a faster calculation

Let's make this a little bit more visual. As an example I'll choose $\ce{(CO2)2}$ and I am using the ridiculous low level of theory HF/3-21G.

First up optimise the equilibrium geometry.

%chk=opt.chk
#P HF/3-21G opt

title required

0 1
C       -1.345383226      0.000000000     -1.107741226
O       -0.793543854     -0.126832008     -0.084830947
O       -1.897222598      0.126832008     -2.130651504
C        1.651522395      0.000000000     -1.243964208
O        2.203361766     -0.126832008     -0.221053930
O        1.099683023      0.126832008     -2.266874487

The result will be the following geometry:

    6
scf done:  -373.128185
 C    -1.546398     0.008833    -1.171593
 O    -0.884661    -0.109039    -0.227710
 O    -2.221473     0.124768    -2.098858
 C     1.852537    -0.008833    -1.180113
 O     2.527612    -0.124768    -0.252847
 O     1.190800     0.109039    -2.123995

Now we want to move the two molecules closer together. If we would use the approach you have come up with, this would be really tedious. First we need to transform the coordinates in a way that we can understand which molecule to move where. Luckily this is not necessary at all. We can simply scan over the "bond length" of the carbon carbon distance. I chose five steps with an increment of -0.2 angstrom.

%chk=modered1.chk
#P HF/3-21G
opt(modredundant)

title

0 1
 C    -1.546398     0.008833    -1.171593
 O    -0.884661    -0.109039    -0.227710
 O    -2.221473     0.124768    -2.098858
 C     1.852537    -0.008833    -1.180113
 O     2.527612    -0.124768    -0.252847
 O     1.190800     0.109039    -2.123995

B 1 4 S 5 -0.2

That will give us the following graph and animation:

You can see that the OCO bond angles changes. The CO bond lengths changes, too, but that is not as visible. This is expected behaviour in a relaxed scan. Only one coordinate is fixed, which is the distance between the carbons. If you are interested in keeping at least the bond angles constant, freeze them:

%chk=modered2.chk
#P HF/3-21G
opt(modredundant)

title

0 1
 C    -1.546398     0.008833    -1.171593
 O    -0.884661    -0.109039    -0.227710
 O    -2.221473     0.124768    -2.098858
 C     1.852537    -0.008833    -1.180113
 O     2.527612    -0.124768    -0.252847
 O     1.190800     0.109039    -2.123995

B 1 4 S 5 -0.2
A 2 1 3 F
A 5 4 6 F

With the result, we'll see again something unexpected happening. There is a sudden rotation at 2.8 angstrom. For the given set of parameters, this is the lowest energy. This is a side effect of a relaxed scan. not everything will run smoothly and sometimes the calculation turns rogue and fails. It really is a delicate matter.

Just for shits and giggles, let's do one that will make us cringe. We will scan also both of the oxygen oxygen distances. Theoretically this should give us 20 data points. Spoiler alert: This will fail at step 8. Try it yourself:

%chk=modered3.chk
#P HF/3-21G
opt(modredundant)

title

0 1
 C    -1.546398     0.008833    -1.171593
 O    -0.884661    -0.109039    -0.227710
 O    -2.221473     0.124768    -2.098858
 C     1.852537    -0.008833    -1.180113
 O     2.527612    -0.124768    -0.252847
 O     1.190800     0.109039    -2.123995

B 1 4 S 5 -0.2
B 2 5 S 5 -0.2
B 3 6 S 5 -0.2

You can see, that with these small molecules it is already hard to predict what happens when you run the calculation. There are traps and pitfalls lining your way. More often than not it is trial and error.

Thankfully a full relaxed scan is not necessary in most cases. Oftentimes a nice little rigid scan will give you a decent result. Have a look at the scan keyword. The only downside to this approach is, that you need a z-matrix to make it work. And putting your brain to that task is also a bit tedious.

I have tried to convert the specifications to a z-matrix putting constants where needed and introducing the scan variable. I'm not sure I succeeded in the way I wanted to, but here we go:

%chk=scan.chk
%nproc=2
%mem=8000MB
#P HF/3-21G scan

title required

0 1
 c
xx   1  1.40000
xx    1  scan    2   90.000
 o    1 oc4         3   90.000      2  88.981
 o    1 oc5         3   90.000      2  -90.00
 c    3  1.93591    1   89.994      4 178.983
xx    6  1.40000    3   90.000      1  90.000
 o    6 oc8         7   90.000      3 180.000
 o    6 oc9         7   90.000      3   0.000

scan        2.794000  5  -0.20
oc4         1.152818
oc5         1.158751
oc8         1.152817
oc9         1.158750

The result is two rods, moving towards each other:

I expect that since you have more than 100 atoms, accomplishing a meaningful scan calculation is a very, very tedious task. It might be much simpler creating the points you want to scan by hand with a molecule editor. But without more information about the actual system, sound advice is difficult.

Answer 2

Firstly, the Gaussian format, while annoyingly unlabeled, has a strict ordering which you can find here.

(Aside from your ordering problem, I highly recommend you check out the Link 0 commands in the Gaussian manual, particularly %OldChk=file which means that it copies the checkpoint file before using it.)

Secondly, I would use (and, actually, have used with an SWNT+arene system) an IRC scan for this: I know it's not quite what you asked about, but if you can start from a reasonably aligned geometry that you think is close to somewhere on the potential valley, it should automatically generate points backward and forward on that valley, without you having to manually specify any coordinates. Using IRC(Report=Read) will spit out the coordinates along with the energy, and from one or two of those sets of coordinates that look close to the minimum, you can perform an optimisation to get the actual minimum. Depending on what you're doing, you may also want to fix many of the nanotube carbons in place before doing this -- it'll dramatically speed up your gradient calculations if I recall correctly. (This is much easier to do if you skip the Z-matrices and just use Cartesian coordinates.)