# Gaussian Convergence Issue for N2+N PES

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?

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

• How did you make you plots? I mean which program do you use? It certainly looks very nice. Good job anyway. ;) – Martin - マーチン Apr 18 '15 at 7:14
• I see a lot of jaggedness in the results. Is this correct? Or is it because of HF method and with CASSCF, things will be much smooth. The reason I have complains about the smoothness is because I would want an analytical fit of the entire PES as the next step. – Tapan Mankodi Jun 2 '15 at 4:46
• @TapanMankodi, we know a priori that HF fails to correctly describe the PES in question, but I'm not sure that the lack of smoothness is the resulf ot the HF failures. I think smoothness also depends on the step used for PES scan which was set to relatively big value of 0.1 A. I would not expect a smooth PES with such settings. – Wildcat Jun 2 '15 at 7:20
• @Wildcat Hmmmm fair enough reasoning about the step size being 0.1A. What recommendation do you have regarding better method should I use. I know that is bit vague question. In other words which should I use that will give me more reliable answers. Also how good cc-pVDZ orbitals? Are they more recommended than STO? Being a mechanical engineering venturing into computational chemistry, your answers are being a great help! – Tapan Mankodi Jun 3 '15 at 6:04
• @TapanMankodi, I'm not a big specialist in the area, but it seems like the system in question requires very advanced electronic structure methods (CASPT2/NEVPT2/MRCI) for its qualitatively correct description. These methods are not even available in Gaussian. I mean, this system is for real experts, not for beginners, you'd better start your venture from something simple. – Wildcat Jun 3 '15 at 9:13

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.

• Yes, Martin, you are perfectly right. HF is hardly appropriate for the calculation in question and spin contamination is indeed quite big. For instance, at a point on PES with $R_1=6.6, R_2=3.7$, SCF converges to a state with $S=1.7416$ which is far away from a quartet state. – Wildcat Apr 17 '15 at 11:37
• So, one should definitely account for electron (Coulomb) correlation, as usual, and should not use HF results by themselves. Doing CASSCF would probably be too expansive, so thinking of the HF calculations as merely an educational exercise, one could do then MP2 or DFT to take electron correlation into account. – Wildcat Apr 17 '15 at 11:41
• @Wildcat To be fair to the method, you are looking for a doublet wave function (when you use the provided input). That, of course, does not make the spin contamination any better. I would not recommend MP2, since you are dealing with nearly degenerate states and the correction to the energy might simply go through the roof, hardly making anything better at that point, maybe even worse. As an educational exercise HF might work, but from the question it is not clear, what the intent behind the study is in the first place. – Martin - マーチン Apr 17 '15 at 11:52
• Sure, sure! To make it more clear: multi-configurational nature of the system in question at long distances can, of course, be dealt with multi-configurational methods only. The meaning of my previous comment was that even at short distances, where 1 CSF is presumably enough, you would better avoid HF and account for Coulomb correlation with MP2/DFT. – Wildcat Apr 17 '15 at 14:49
• @Tapan Basically yes. In the line above the geometry you indicate charge (first number) and multiplicity (second number) as 2S+1. – Martin - マーチン Jun 2 '15 at 6:35