Issues with gamess external basis set d2tzvp and iodine

[rpseng]

Hello, thanks a lot for this repository and the instructions on how to use EXTBASIS.txt with gamess.

All works fine with d2tzvp except for molecules with iodine, e.g. iodomethane. The energy does not converge, keeps oscilating:

ITER EX DEM     TOTAL ENERGY        E CHANGE  DENSITY CHANGE     ORB. GRAD
   1  0  0    -2688.9031574371 -2688.9031574371  47.972062755   0.000000000
   2  1  0    -2712.8674452076   -23.9642877705 166.206976187   0.869817704
   3  2  0    -2505.6269175992   207.2405276084 165.606852333   1.502258158
   4  3  0    -2670.6303385744  -165.0034209752  65.987403521   0.866474212
   5  0  0    -2466.3695781028   204.2607604716  28.721711243   2.232559846
   6  1  0    -2713.5127307253  -247.1431526224 159.015100570   0.971159474
   7  2  0    -2505.1464430413   208.3662876839 157.638077469   1.579522623
   8  3  0    -2670.5450082936  -165.3985652522  65.781634675   0.879239380
   9  0  0    -2466.3019640269   204.2430442667  28.730265657   2.234209766
  10  1  0    -2713.5157210149  -247.2137569879 159.102230658   0.972023360
  11  2  0    -2505.1475679112   208.3681531036 157.727222664   1.579171900
  12  3  0    -2670.5457304377  -165.3981625265  65.783387844   0.879128294
  13  0  0    -2466.3025277625   204.2432026752  28.730238381   2.234197719
  14  1  0    -2713.5156762485  -247.2131484861 159.101489561   0.972018067
  15  2  0    -2505.1475500997   208.3681261488 157.726464386   1.579174814
  16  3  0    -2670.5457235111  -165.3981734114  65.783371216   0.879129344
  17  0  0    -2466.3025223706   204.2432011406  28.730238749   2.234197835
  18  1  0    -2713.5156766288  -247.2131542582 159.101496645   0.972018122
  19  2  0    -2505.1475502493   208.3681263795 157.726471635   1.579174786
  20  3  0    -2670.5457235752  -165.3981733259  65.783371371   0.879129334
  21  0  0    -2466.3025224205   204.2432011547  28.730238745   2.234197834
  22  1  0    -2713.5156766248  -247.2131542043 159.101496581   0.972018121
  23  2  0    -2505.1475502481   208.3681263767 157.726471569   1.579174786
  24  3  0    -2670.5457235745  -165.3981733264  65.783371370   0.879129334
  25  0  0    -2466.3025224199   204.2432011546  28.730238745   2.234197834
  26  1  0    -2713.5156766249  -247.2131542049 159.101496570   0.972018121
  27  2  0    -2505.1475502480   208.3681263769 157.726471558   1.579174786
  28  3  0    -2670.5457235745  -165.3981733265  65.783371369   0.879129334
  29  0  0    -2466.3025224200   204.2432011545  28.730238745   2.234197834
  30  1  0    -2713.5156766247  -247.2131542047 159.101496566   0.972018121

When checking the basis set exchange website for this atom, there is an extra ECP section:
https://www.basissetexchange.org/basis/def2-tzvp/format/gamess_us/?version=1&elements=53

I'm not sure if this is the source of the problem or if this ECP section should be included somehow. A minimal input file follows below:

 $contrl scftyp=rhf $end
 $contrl mplevl=0 $end

 $system timlim=10000 mwords=20 $end
 $contrl ISPHER=1 $END

! https://github.com/shoubhikraj/molecular-modelling/tree/main/basis-sets
 $BASIS GBASIS=d2tzvp EXTFIL=.t. $END

 $guess guess=HUCKEL $end
 $SYSTEM MWORDS=500 $END

 $DATA
http://cactus.nci.nih.gov/chemical/structure/IODOMETHANE/file?format=sdf&get3d=True
C1
I     53.0     -0.2367000000    0.0000000000    0.0000000000
C      6.0      1.9253000000    0.0000000000    0.0000000000
H      1.0      2.2886000000    0.0000000000    1.0277000000
H      1.0      2.2886000000   -0.8900000000   -0.5138000000
H      1.0      2.2886000000    0.8900000000   -0.5138000000
 $END

Any insight on this would be very much appreciated.

[shoubhikraj]

@rpseng Hi, glad to hear the files are useful 😄.

You guessed correctly that the problem is coming from the ECP. From what I can see in basis set exchange, all elements after Kr have ECP's. So it's not surprising that the calculation fails. Adding the $ECP group for iodine should fix your problem:

$contrl scftyp=rhf $end
 $contrl mplevl=0 $end

 $system timlim=10000 mwords=20 $end
 $contrl ISPHER=1 PP=READ $END

! https://github.com/shoubhikraj/molecular-modelling/tree/main/basis-sets
 $BASIS GBASIS=d2tzvp EXTFIL=.t. $END

 $guess guess=HUCKEL $end
 $SYSTEM MWORDS=500 $END

 $DATA
http://cactus.nci.nih.gov/chemical/structure/IODOMETHANE/file?format=sdf&get3d=True
C1
I     53.0     -0.2367000000    0.0000000000    0.0000000000
C      6.0      1.9253000000    0.0000000000    0.0000000000
H      1.0      2.2886000000    0.0000000000    1.0277000000
H      1.0      2.2886000000   -0.8900000000   -0.5138000000
H      1.0      2.2886000000    0.8900000000   -0.5138000000
 $END

 $ECP
I-ECP GEN    28    3
4     ----- f-ul potential -----
    -21.84204000      2      19.45860900
    -28.46819100      2      19.34926000
     -0.24371300      2       4.82376700
     -0.32080400      2       4.88431500
7     ----- s-f potential -----
     49.99429300      2      40.01583500
    281.02531700      2      17.42974700
     61.57332600      2       9.00548400
     21.84204000      2      19.45860900
     28.46819100      2      19.34926000
      0.24371300      2       4.82376700
      0.32080400      2       4.88431500
8     ----- p-f potential -----
     67.44284100      2      15.35546600
    134.88113700      2      14.97183300
     14.67505100      2       8.96016400
     29.37566600      2       8.25909600
     21.84204000      2      19.45860900
     28.46819100      2      19.34926000
      0.24371300      2       4.82376700
      0.32080400      2       4.88431500
10    ----- d-f potential -----
     35.43952900      2      15.06890800
     53.17605700      2      14.55532200
      9.06719500      2       6.71864700
     13.20693700      2       6.45639300
      0.08933500      2       1.19177900
      0.05238000      2       1.29115700
     21.84204000      2      19.45860900
     28.46819100      2      19.34926000
      0.24371300      2       4.82376700
      0.32080400      2       4.88431500
C-ECP NONE
H-ECP NONE
H-ECP NONE
H-ECP NONE
 $END

GAMESS's ECP input is a bit tricky in that you need to define the ECP for each atom in the exact same order they appear in $DATA, and you need to do it for all the atoms - even the ones that do not have ECP. (https://www.msg.chem.iastate.edu/gamess/GAMESS_Manual/docs-input.txt). Once you put the ECP in for Iodine, the calculation runs smoothly.

         --------------------------
                RHF SCF CALCULATION
         --------------------------

    NUCLEAR ENERGY =        60.9020337474
    MAXIT =   30     NPUNCH=    2
    EXTRAP=T  DAMP=F  SHIFT=F  RSTRCT=F  DIIS=F  DEM=F  SOSCF=T
    DENSITY MATRIX CONV=  1.00E-05
    SOSCF WILL OPTIMIZE    1394 ORBITAL ROTATIONS, SOGTOL=   0.250
    MEMORY REQUIRED FOR RHF ITERS=    126541 WORDS.

ITER EX DEM     TOTAL ENERGY        E CHANGE  DENSITY CHANGE     ORB. GRAD
  1  0  0     -325.7821902639  -325.7821902639   2.334790141   0.000000000
         ---------------START SECOND ORDER SCF---------------
  2  1  0     -336.1580140330   -10.3758237691   0.450395631   0.155375000
  3  2  0     -336.2542577709    -0.0962437380   0.126198268   0.058815788
  4  3  0     -336.2870328290    -0.0327750581   0.020510540   0.007031238
  5  4  0     -336.2875663758    -0.0005335468   0.008046922   0.003710974
  6  5  0     -336.2877495081    -0.0001831323   0.001741688   0.000532772
  7  6  0     -336.2877606451    -0.0000111370   0.000878706   0.000358150
  8  7  0     -336.2877619997    -0.0000013546   0.000318105   0.000091646
  9  8  0     -336.2877621418    -0.0000001421   0.000064803   0.000020826
 10  9  0     -336.2877621527    -0.0000000109   0.000011581   0.000005576
 11 10  0     -336.2877621532    -0.0000000005   0.000003352   0.000001168
 12 11  0     -336.2877621532    -0.0000000000   0.000000912   0.000000215

         -----------------
         DENSITY CONVERGED
         -----------------
    TIME TO FORM FOCK OPERATORS=       0.5 SECONDS (       0.0 SEC/ITER)
    TIME TO SOLVE SCF EQUATIONS=       0.0 SECONDS (       0.0 SEC/ITER)

FINAL RHF ENERGY IS     -336.2877621532 AFTER  12 ITERATIONS

[rpseng]

Hello @shoubhikraj, thanks a lot for the prompt response!

Now I'm thinking on how I can automate this. I'm planning to update the basis set we use on this project https://github.com/lvpp/sigma, but we need to run it for about 2000 molecules. Current automation uses openbabel to create the gamess input from a .mol and an extra file with the gamess keys we use. With this per atom extra input, that will be very tricky...

[shoubhikraj]

@rpseng Yes automation would be tricky, but doable I believe. You are using python for scripting, from what I can see? Then you could read the generated input and get all atoms in order, and then have another text file with all the ECP's for different atoms, and put the proper ECP's together to write the $ECP block in input file.

Unfortunately, GAMESS has not put def2 bases in its library and the input style for ECP is awkward, so without this kind of python based hack I don't think there is a good solution.