Gaussian Convergence Issue for N2+N PES

Question

I am trying to obtain a potential energy surface (PES) for $\ce{N2 + N}$ combination using scan feature of Gaussian. I have tried different combinations of basis sets and methods but unfortunately at some point in the iteration, the convergence criteria is not met and I get an "Aborted Core dump" closure in the log file. The Gaussian code is:

# HF/6-31G test scan 

NNN PES

0 2
N
N 1 R1
N 2 R2 1 120.

R1 0.5 80 0.1
R2 0.5 80 0.1

Where am I going wrong. More importantly, I would be really thankful if someone could suggest how to choose a method and basis set for a type of problem. I was thinking if I can get some answer using HF/6-31G, I could use this as an input for a higher order method like cassf. What that be a good idea?

Answer 1

You did not specify at which geometry you have convergence problems, so at first I thought it happens when all three $\ce{N}$ atoms are far away from each other. But then I run the calculation by myself and found out that it is not the case: the SCF convergence failure happens at R1=1.2 and R2=1.1. At this geometry SCF starts to oscillate after 20 iterations back and force as shown in a graph below, so simply increasing the number of SCF iteration would not help.

For that reason I quickly checked few other well known tricks to avoid the problems with SCF convergence (see, for instance, here) and I recommend to use a quadratically convergent SCF (QC-SCF) by adding scf=qc to the route section:

# HF/6-31G test scan scf=qc

Note though that while QC-SCF is more reliable than the default DIIS algorithm it is also slower, so expect an increase in the computation time. For that reason, there is also an interesting option scf=xqc which first tries the default algorithm and switches to a quadratically convergent one only if the default has not converged.

# HF/6-31G test scan scf=xqc

---

However, switching to a quadratically convergent algorithm did not help, since SCF still fails to converge but at some other geometries. So, I investigated the matter more closely, and in addition to scf=xqc I also recommend to improve both radial and angular flexibility of you basis set by adding diffuse and polarization functions into it

# HF/6-31+G(d,p) test scan scf=xqc

---

After running the calculation for an hour it reached 2656 points (which is about 40% of all points) without SCF convergence problems, but in case if this variant still fails:

• Consider lowering down the convergence criteria by changing scf=xqc to scf=(xqc,conver=n) with n being, say, 6 (it is 8 by default).
• If failures happen when at least one $\ce{N}$ is far away from two others, consider lowering down the maximum distance for your PES scan. In my opinion 8.5Å is way to big, 5Å would probably be more than enough.

---

I finally succeeded with the following route section

# HF/6-31+G(d,p) Scan SCF=(XQC,Conver=6)

Here are the results and below is the plot of the PES. :D

Answer 2

Apart from the technical aspects, that Wildcat already provided, I would assume, that the problem is far too complicated to be treated with Hartree-Fock correctly. It might be even not possible.

Hartree-Fock is a method that describes no coulomb correlation at all. I would assume, that this is vital to the system in any case.

At long distances you will have the problem, that you cannot describe the system in a doublet state any more. Nitrogen as an atom has five valence electrons, and according to Hund's rule it should be a quartet spin state. So if you have the separation $\ce{N + N2}$ this should also be a quartet state. In this case you will have (nearly) degenerate orbitals, that will cause the SCF to oscillate.

When you come to the situation $\ce{N + N + N}$, you are effectively treating it as a super-molecule, where you couple three quartet states. The only way of qualitatively treating this correctly would be a complete active space ansatz. Here again it is most likely, that the SCF will oscillate because of (nearly) degenerate orbitals.

Since you are freezing the bond angle at 120.0°, you should also be aware, that the equilibrium geometry will not be part of your designed potential energy surface. I would assume, that it is linear, analogous to the azide ion.

In any case I would expect the wave function to be highly spin contaminated.

Running a scan at the CASSCF level of theory is very time consuming and involves high effort computations. I would guess you need to pick at least 15 active electrons in 15 orbitals. For starters it might work only treating the π system with CAS, but that still leaves 9 electrons in 6 orbitals.