Atomic start guess

Initial SCF run

We consider the Hartree-Fock calculation of porphyrin with a mercury center

alternate text

The SCF input file HF.inp reads:

**DIRAC
.WAVE FUNCTION
.ANALYZE
**WAVE FUNCTIONS
.SCF
*SCF
.MAXITR
20
.PRINT
2
.CLOSED SHELL
120 120
**ANALYZE
.MULPOP
*END OF

for an SCF calculation with a Mulliken population analysis. We specify the maximum number of iterations allowed for th caclulation to 20 and as we are chosing a molecule with inversion symmetry we specify the occupation number for both irreducible representations under the keyword CLOSED SHELL. The print level is set to 2 in order to have the eigenvlues printed for each iteration. The molecular input hg-porphyrin.mol reads :

BASIS
6-31 and finite DK3 basis sets, geometry in au, ref ECP: Theor.Chim.Acta 1990,77,123 (from Gaussian)
--------------------------------------------------------
C   4    3  X Y Z
        7.    2
N1    0.000000        4.1196801826    0.000000
N2    4.1196801826    0.000000        0.000000
LARGE BASIS 6-31G
        6.    5
C1       4.642960672    4.642960672    0.000000
C2       2.127242869    5.577562894    0.000000
C3       1.2976452113   8.20576858026  0.000000
C4       8.20576858026  1.2976452113   0.000000
C5       5.577562894    2.127242869    0.000000
LARGE BASIS 6-31G
        1.    3
H1        6.1040918256  6.1040918256   0.000000
H2        2.536971561   9.84908074236  0.000000
H3        9.84908074236 2.536971561    0.000000
LARGE BASIS 6-31G
       80.    1
Hg    0.000000     0.000000     0.000000
LARGE    4    1    1    1    1
F  23   0    0                                        # s-functions
   5.4107818D+07
   1.2757417D+07
   3.8200687D+06
   1.2601464D+06
   4.4355086D+05
   1.6477405D+05
   6.4498870D+04
   2.6536057D+04
   1.1391885D+04
   5.0460108D+03
   2.2877229D+03
   1.0547293D+03
   4.9377065D+02
   2.3032535D+02
   1.1057622D+02
   4.7454363D+01
   2.4814292D+01
   1.0508682D+01
   5.3690964D+00
   1.7593797D+00
   8.2178073D-01
   1.6991016D-01
   5.9613310D-02
F  23   0    0                                        # p-functions
   2.1202540D+07
   4.1579701D+06
   1.0372139D+06
   2.9906582D+05
   9.5123364D+04
   3.3364166D+04
   1.2858845D+04
   5.3680904D+03
   2.3907706D+03
   1.1241754D+03
   5.5303069D+02
   2.8249864D+02
   1.4871000D+02
   7.9474383D+01
   4.3674317D+01
   2.4527063D+01
   1.4058201D+01
   7.9450158D+00
   4.5067591D+00
   2.4440572D+00
   1.2660103D+00
   6.1522560D-01
   2.5113161D-01
F  15   0    0                                        # d-functions
   1.4315825D+04
   3.7250818D+03
   1.3315379D+03
   5.5776192D+02
   2.5656959D+02
   1.2554901D+02
   6.3698254D+01
   3.3232923D+01
   1.7367497D+01
   8.9182853D+00
   4.5159858D+00
   2.1444878D+00
   9.7876159D-01
   4.1519292D-01
   1.5886171D-01
F  10   0    0                                        # f-functions
   1.3065007D+03
   4.4185000D+02
   1.8456705D+02
   8.5884633D+01
   4.2090822D+01
   2.1327466D+01
   1.0817659D+01
   5.3881551D+00
   2.5585604D+00
   1.0775039D+00
FINISH

The calculations are carried out in D2h symmetry, so that the Fock matrix is real. From the output file we read:

Atoms and basis sets
--------------------

Number of atom types:     4
Total number of atoms:   37

  label    atoms   charge   prim    cont     basis
  ----------------------------------------------------------------------
  N2          4       7      22       9      L  - [10s4p|3s2p]
                             58      58      S  - [4s10p4d|4s10p4d]
  C5         20       5      22       9      L  - [10s4p|3s2p]
                             58      58      S  - [4s10p4d|4s10p4d]
  H3         12       1       4       2      L  - [4s|2s]
                             12      12      S  - [4p|4p]
  Hg          1      80     282     282      L  - [23s23p15d10f|23s23p15d10f]
                            635     635      S  - [23s38p33d15f10g|23s38p33d15f10g]
  ----------------------------------------------------------------------
                            858     522      L  - large components
                           2171    2171      S  - small components
  ----------------------------------------------------------------------
  total:     37     220    3029    2693

showing that in AO-basis the Fock-matrix has the dimension 2693 x 2693 and thus requires 7252249 words. The SCF run will require something like three Fock matrices, so in this case 21 Mwords.

Unfortunately, this HF does not converge:

                                   SCF - CYCLE
                                   -----------

--------------------------------------------------------------------------------------------------------------------------------
           Energy               ERGVAL    FCKVAL    EVCVAL      Conv.acc    CPU          Integrals   Time stamp
--------------------------------------------------------------------------------------------------------------------------------
It.    1    -7794.747305914      0.00D+00  0.00D+00  0.00D+00              10.89634323s   Scr. nuclei    Fri Oct 19
It.    2    -20119.30505450      1.23D+04 -9.63D+01  1.29D+02               4min 6.910s   LL SL          Fri Oct 19
It.    3    -17230.23561991     -2.89D+03  1.95D+02  2.69D+02   DIIS   2    3min55.989s   LL SL          Fri Oct 19
It.    4    -20286.13741996      3.06D+03 -8.61D+01  9.85D+01   DIIS   3    3min42.372s   LL SL          Fri Oct 19
It.    5    -20126.54429484     -1.60D+02 -9.61D+01  4.27D+01   DIIS   4    3min52.176s   LL SL          Fri Oct 19
It.    6    -18673.18452215     -1.45D+03  1.36D+02  1.55D+02   DIIS   5    3min48.236s   LL SL          Fri Oct 19
It.    7    -18727.77909394      5.46D+01 -1.82D+00  1.48D+02   DIIS   6    3min53.450s   LL SL          Fri Oct 19
It.    8    -18755.21459721      2.74D+01 -8.49D-01  1.45D+02   DIIS   7    4min 4.378s   LL SL          Fri Oct 19
It.    9    -18760.27704618      5.06D+00 -1.58D-01  1.45D+02   DIIS   8    3min40.865s   LL SL          Fri Oct 19
It.   10    -18773.33929407      1.31D+01 -4.09D-01  1.44D+02   DIIS   9    3min51.322s   LL SL          Fri Oct 19
It.   11    -18779.70838778      6.37D+00 -1.28D-01  1.44D+02   DIIS   9    3min42.698s   LL SL          Fri Oct 19
It.   12    -18953.17392936      1.73D+02 -4.17D+00  1.33D+02   DIIS   9    3min52.510s   LL SL          Fri Oct 19
It.   13    -18938.57698783     -1.46D+01  5.29D-01  1.34D+02   DIIS   9    3min43.147s   LL SL          Fri Oct 19
It.   14    -18744.84737430     -1.94D+02  8.55D+00  1.59D+02   DIIS   9    3min56.082s   LL SL          Fri Oct 19
It.   15    -18164.22290167     -5.81D+02 -2.01D+02  4.93D+01   DIIS   9    3min48.373s   LL SL          Fri Oct 19
It.   16    -16110.21649469     -2.05D+03  2.64D+02  2.90D+02   DIIS   9    3min47.714s   LL SL          Fri Oct 19
It.   17    -16272.74479238      1.63D+02 -4.12D+00  2.75D+02   DIIS   9    3min49.219s   LL SL          Fri Oct 19
It.   18    -16280.64616752      7.90D+00 -2.09D-01  2.74D+02   DIIS   9    3min40.992s   LL SL          Fri Oct 19
It.   19    -16279.46291382     -1.18D+00  5.93D-02  2.75D+02   DIIS   9    3min37.488s   LL SL          Fri Oct 19
It.   20    -16289.37827128      9.92D+00 -2.58D-01  2.74D+02   DIIS   9    3min42.223s   LL SL          Fri Oct 19

We see that the corrected bare nucleus start is quite far from the energies of subsequent iterations. Although the energy in the second iteration drops considerably it bounces back and gets stuck. We will therefore switch to an atomic start [vanLenthe2006] . Observe that this caclulation will take more than an hour.

Preparing the atomic start

For each atom we now run a SCF calculation in full (linear) symmetry based on the atomic ground state configuration. For open-shell atoms this amounts to an average-of-configuration (AOC) calculation at the HF level, and a fractional occupation calculation at the DFT level.

Nitrogen atom

The input file N.inp reads:

**DIRAC
.WAVE FUNCTION
.ANALYZE
**GENERAL
.ACMOUT
**WAVE FUNCTIONS
.SCF
*SCF
.CLOSED SHELL
4 0
.OPEN SHELL
1
3/0,6
**ANALYZE
.MULPOP
*END OF

Note the keywords ACMOUT asking that the MO coeffcients transformed to no symmetry and written to file DFACMO. The open shell input specifies 1 open shell with 3 electrons in 0 and 6 spinors respectively for the two irriducible representations. The molecular input N.mol is

BASIS
6-31 and finite DK3 basis sets, geometry in au, ref ECP: Theor.Chim.Acta 1990,77,123 (from Gaussian)
--------------------------------------------------------
C   1
        7.    1
N      0.0 0.0   0.0
LARGE BASIS 6-31G
FINISH

The minimal pam command is

pam --mol=N.mol --inp=N.inp --get=DFACMO
mv DFACMO ac.N

Other atoms

The input file for the other atoms are prepared in similar manner. For carbon we use the occupation

.CLOSED SHELL
4 0
.OPEN SHELL
1
2/0,6

For mercury we use

.CLOSED SHELL
42 38

whereas the hydrogen atom is run as a one-electron system, meaning that the keyword for ignoring two-particle interactions ONESYS will be invoked

**DIRAC
.WAVE FUNCTION
.ANALYZE
**GENERAL
.ACMOUT
**HAMILTONIAN
.ONESYS
**WAVE FUNCTIONS
.SCF
*SCF
.CLOSED SHELL
2 0
**ANALYZE
.MULPOP
*END OF

(The occupation given above is simply to control orbitals printed in the subsequent Mulliken population analysis).

All atoms can be run in a single shot, e.g. using bash

for atom in H C N Hg
do
  pam --mol=${atom}.mol --inp=${atom}.inp --get=DFACMO
  mv DFACMO ac.${atom}
done

We are now ready to run the atomic start.

Running the atomic start

The input file atomstart.inp reads

**DIRAC
.WAVE FUNCTION
.ANALYZE
**WAVE FUNCTIONS
.SCF
*SCF
.ATOMST
AFNXXX 2
1,2
1.00
3..5
0.5
AFCXXX 2
1,2
1.00
3..5
0.33
AFHXXX 1
1
0.5
AFHGXX 1
1..40
1.00
.CLOSED SHELL
 120 120
**ANALYZE
.MULPOP
*END OF

The keyword ATOMST is followed by input for each atomic type. The first input AFYYXX specifies the file to be read for an atomtype YY, and next to it a number specifying inversion symmtry (2) or not (1). Next, for each symmetry, the spinors of that symmetry is specified followed by the occupation of the spinors. The occupations chosen corresponds to those of the atomic runs, but the user may modify this at will. Please note that the order of atoms corresponds to the order they appear in the molecule file.

Before running the calculation the user must make provide links to the atomic coefficient files

ln -s ac.N AFNXXX
ln -s ac.C AFCXXX
ln -s ac.H AFHXXX
ln -s ac.Hg AFHGXX

and the calculation is then run as

pam --mol=hg-porphyrin.mol --inp=atomstart.inp --copy="AF*" --outcmo

This calculation converges smoothly after 17 iterations

                                   SCF - CYCLE
                                   -----------

* Convergence on norm of error vector (gradient).
  Desired convergence:1.000D-07
  Allowed convergence:1.000D-06

* ERGVAL - convergence in total energy
* FCKVAL - convergence in maximum change in total Fock matrix
* EVCVAL - convergence in error vector (gradient)
--------------------------------------------------------------------------------------------------------------------------------
           Energy               ERGVAL    FCKVAL    EVCVAL      Conv.acc    CPU          Integrals   Time stamp
--------------------------------------------------------------------------------------------------------------------------------
It.    1    -20645.07067749      2.06D+04  0.00D+00  0.00D+00   Atomic s   32.88100123s   LL SL          Mon Mar  4
It.    2    -20633.15458108     -1.19D+01  6.07D-01  8.19D-01               1min11.442s   LL SL          Mon Mar  4
It.    3    -20633.38854493      2.34D-01 -2.13D-01  3.24D-01   DIIS   2    1min 9.861s   LL SL          Mon Mar  4
It.    4    -20633.41415475      2.56D-02  1.52D-01  1.72D-01   DIIS   3    1min 3.733s   LL SL          Mon Mar  4
It.    5    -20633.42153720      7.38D-03 -6.24D-02  2.79D-02   DIIS   4    1min 2.223s   LL SL          Mon Mar  4
It.    6    -20633.42175881      2.22D-04 -4.58D-03  8.00D-03   DIIS   5    1min 3.556s   LL SL          Mon Mar  4
It.    7    -20633.42178693      2.81D-05  2.31D-03  3.88D-03   DIIS   6    1min 1.382s   LL SL          Mon Mar  4
It.    8    -20633.42179098      4.05D-06 -1.68D-03  1.75D-03   DIIS   7   57.84219360s   LL SL          Mon Mar  4
It.    9    -20633.42179253      1.55D-06  3.91D-04  5.46D-04   DIIS   8   57.79122925s   LL SL          Mon Mar  4
It.   10    -20633.42179283      3.01D-07 -5.63D-05  2.16D-04   DIIS   9   55.75250244s   LL SL          Mon Mar  4
It.   11    -20633.42179286      3.30D-08  1.89D-05  7.63D-05   DIIS   9   53.37689209s   LL SL          Mon Mar  4
It.   12    -20633.42179287      4.12D-09 -1.74D-05  2.92D-05   DIIS   9   51.99108887s   LL SL          Mon Mar  4
It.   13    -20633.42179287      7.42D-10  6.49D-06  8.29D-06   DIIS   9   51.30120850s   LL SL          Mon Mar  4
It.   14    -20633.42179287      3.13D-10 -1.22D-06  2.61D-06   DIIS   9   48.06268311s   LL SL          Mon Mar  4
It.   15    -20633.42179287     -2.18D-10  5.07D-07  5.08D-07   DIIS   9   48.58166504s   LL SL          Mon Mar  4
It.   16    -20633.42179287     -2.15D-10 -2.13D-07  2.33D-07   DIIS   9   43.13543701s   LL SL          Mon Mar  4
It.   17    -20633.42179287      5.42D-10  1.06D-07  9.96D-08   DIIS   9   42.37854004s   LL SL          Mon Mar  4
--------------------------------------------------------------------------------------------------------------------------------
* Convergence after   17 iterations.