# Why does the SCF not converge in a structure optimization with GAMESS of the thiourea molecule at the PM6 level of theory?

I'm using GAMESS to run a geometry optimization on a thiourea molecule, using PM6. However, the output file keeps displaying the error message, "NO FORCE FIELD, SCF DOES NOT CONVERGE." Is there an error in my input code?

 $$CONTRL SCFTYP=RHF RUNTYP=OPTIMIZE MAXIT=100 MULT=1 COORD=UNIQUE NOSYM=1$$END
$$BASIS GBASIS=PM6$$END
$$DATA THIOUREA C1 C 6 -0.000001 -0.216235 -0.000006 N 7 -1.149850 -0.964256 -0.062205 N 7 1.149844 -0.964254 0.062165 S 16 0.000003 1.474558 -0.000003 C 6 -2.468377 -0.327323 0.016251 C 6 2.468378 -0.327328 -0.016221 H 1 -3.255938 -1.092284 -0.077195 H 1 -2.599975 0.405207 -0.806499 H 1 -2.609097 0.216356 0.970441 H 1 3.255932 -1.092302 0.077168 H 1 2.599974 0.405143 0.806581 H 1 2.609108 0.216419 -0.970372 H 1 -1.133325 -1.963288 0.067490 H 1 1.133319 -1.963297 -0.067428$$END

• In these cases, the first thing to do might be to increase the number of SCF cycles, so try MAXIT=150, and see what happens.
– Pier
Jul 12 '19 at 0:22
• I am unfamiliar with GAMESS, but make sure you have specified the method correctly. Also, check your input geometry visually and consider switching to another method. Maybe the C-S bond is badly described in PM6. Jul 12 '19 at 5:48
• I couldn't find any indication that pm6 is actually available in GAMESS. Maybe try pm3, and see whether it works. I think pm6 would do fine for the molecule, pm3 probably too. If you require pm6 you might have to try a different ESS code. Jul 12 '19 at 11:28
• I spent some time yesterday trying it on WebMo (I have no local GAMESS installations on the OS I am working on right now), and even if I am not familiar with GAMESS, it looks like for some weird reason PM6 doesn't like that molecule (increased optimization and scf cycles, modified convergence criteria). If I recall correctly, anyway, it worked with PM3. Jul 13 '19 at 20:31
• I believe the problem has to lie within GAMESS. I just ran an optimisation with Gaussian 16 B.01 and pm6, which looks like you'd would expect it. ( cc @The_Vinz ) Could you please specify the version you are using. Have you tried other molecules; did they work? Maybe good choices to localise the problem would be urea, CS2/CS, maybe the problem lies there. I'll see if I can get GAMESS to work, then I might try again. Jul 23 '19 at 11:47

TL;DR: PM6 only works for elements up to fluoride in GAMESS.

I have played around with the GAMESS version installed at our cluster, the ancient 5 DEC 2014 (R1) and I came to the conclusion that there must be a 'bug' in the implementation of PM6 with GAMESS. (I reproduced the problem with 20 APR 2017 (R1); I have also downloaded 30 Jun 2019 (R1), but I won't try it as I think the result is the same.) When I tried it with Gaussian 16 Rev. B.01 on the PM6 level of theory it converged fine.

Actually, I have looked at the documentation, and it seems to be officially not supported (from INTRO.DOC):

Capabilities

[...]
3. Calculates semi-empirical MNDO, AM1, or PM3 models
using RHF, UHF, ROHF, or GVB wavefunctions.
[...]


Or from INPUT.DOC:

        * * * semiempirical basis sets * * *

GBASIS = MNDO - selects MNDO model Hamiltonian
= AM1  - selects AM1 model Hamiltonian
= PM3  - selects PM3 model Hamiltonian
= RM1  - selects RM1 model Hamiltonian
= DFTB - selects tight binding Hamiltonian


I am pretty sure the error is associated with sulfur. I have tried a few molecules, among them $$\ce{H3C-NH2}$$, $$\ce{(H3C)2NH}$$, $$\ce{(H3C)2CO}$$ (see input below), urea, which converged fine. Whenever I tried a molecule containing sulfur the SCF seemed to be oscillating. I have tested $$\ce{CS2}$$, $$\ce{SO2}$$, $$\ce{(H3C)2CS}$$, and (obviously) thiourea. (With PM3 everything worked fine.)

My best guess is that they simply didn't finish implementing it. You could issue a support ticket on GitHub: issues. I found an article by (the here known) Jan Jensen from 2013: Project status: PM6 and PCM in GAMESS. And then from 2014: PM6, DFTB, GAMESS, enzyme design, protein-ligand binding, and protein structure validation. I stopped there.

I'd recommend using a different program, the obvious choice would be MOPAC (http://openmopac.net/). Gaussian has it implemented and a couple others probably too. Or you switch to a completely different method like xtb, HF-3c (in ORCA), DFTB (e.g. DFTB+), etc.. Depending on your system size, you can probably also get DFT based calculations done.

Here is also something related: GAMESS crash course?

 $$CONTRL SCFTYP=RHF RUNTYP=OPTIMIZE MAXIT=100 MULT=1 COORD=UNIQUE NOSYM=1$$END
$$BASIS GBASIS=PM6$$END
$$DATA C3H6O C1 CARBON 6.0 0.156332162 -0.059231571 0.417752882 OXYGEN 8.0 1.200569512 -0.535604279 0.831319560 CARBON 6.0 -0.037914597 1.466530485 0.340822109 HYDROGEN 1.0 0.279389785 1.860652999 -0.631337255 HYDROGEN 1.0 -1.089141584 1.742544046 0.481697509 HYDROGEN 1.0 0.545018653 1.982387227 1.112307247 CARBON 6.0 -1.053925162 -0.882111362 -0.061565447 HYDROGEN 1.0 -0.987243509 -1.102851552 -1.132983096 HYDROGEN 1.0 -1.121341433 -1.838731138 0.468944868 HYDROGEN 1.0 -1.993938204 -0.344384630 0.106761952$$END