4.3 Generating MO Start Vectors

4.3.1 The MO Start Vectors Menu

This menu serves to define the occupation numbers, and to generate the start vectors for HF-SCF and DFT calculations. They may be constructed from earlier SCF calculations (perhaps employing another basis set, type use), by Hamilton core guess (hcore), or by an extended Hückel calculation which can be performed automatically (eht). An occupation of the start orbitals will be proposed and can be modified if desired.

OCCUPATION NUMBER & MOLECULAR ORBITAL DEFINITION MENU  
 
 CHOOSE COMMAND  
 infsao     : OUTPUT SAO INFORMATION  
 atb        : Switch for writing MOs in ASCII or binary format  
 eht        : PROVIDE MOS && OCCUPATION NUMBERS FROM EXTENDED HUECKEL GUESS  
 use <file> : SUPPLY MO INFORMATION USING DATA FROM <file>  
 man        : MANUAL SPECIFICATION OF OCCUPATION NUMBERS  
 hcore      : HAMILTON CORE GUESS FOR MOS  
 flip       : FLIP SPIN OF A SELECTED ATOM  
 &          : MOVE BACK TO THE ATOMIC ATTRIBUTES MENU  
 THE COMMANDS  use  OR  eht  OR  *  OR q(uit) TERMINATE THIS MENU !!!  
 FOR EXPLANATIONS APPEND A QUESTION MARK (?) TO ANY COMMAND

Recommendation

You will normally only need to enter eht. For the EHT-guess, define will ask for some specifications, and you should always choose the default, i.e. just <enter>. Of importance is only the molecular charge, specified as 0 (neutral, default), 1 or -1 etc.

Based on the EHT orbital energies define proposes an occupation. If you accept you are done, if not you get the “occupation number assignment menu” explained in 4.3.2.

Description of Commands

infsao

Command infsao provides information about the symmetry adapted basis which is used for the SCF-calculation. To exploit the molecular symmetry as efficiently as possible, TURBOMOLE programs do not use the simple basis which you specified during your define session. Instead it builds linear combinations of equal basis functions on different—but symmetry equivalent—atoms. This basis is then called the SAO-basis (Symmetry Adapted Orbital). It has the useful property that each basis function transformed to this basis transforms belongs to one irreducible representation of the molecular point group (that is, the basis reflects the full molecular symmetry as specified by the Schönflies symbol). infsao gives you a listing of all symmetry adapted basis functions and their constituents either on file or on the screen. This may help you if you want to have a closer look at the SCF vectors, because the vector which is output by program dscf is written in terms of these SAOs.

atb

Molecular orbitals can be written either in ASCII or in binary format. This command switches from one option to the other, and it is highly recommended to read which setting is currently active. ASCII format is portable and allows the usage of TURBOMOLEinput files on different systems with incompatible binary format. Binary format is faster and smaller files will be written. The external program atbandbta can be used to transform existing mos, alpha, and beta files from ASCII to binary format and vice versa.

eht

eht performs an extended Hückel calculation for your molecule. The orbital energies available from this calculation are then used to provide occupation numbers for your calculation and the Hückel MOs will be projected onto the space that is spanned by your basis set. This start-vectors are not as good as the ones which may be obtained by projection of an old SCF vector, but they are still better than the core Hamiltonian guess that is used if no start vectors are available. When using this command, you will be asked if you want to accept the standard Hückel parameters and to enter the molecular charge. Afterwards you will normally get a list of the few highest occupied and lowest unoccupied MOs, their energies and their default occupation. If you don’t want to accept the default occupation you will enter the occupation number assignment menu, which is described in Section 4.3.2. Note that the occupation based on the Hückel calculation may be unreliable if the difference of the energies of the HOMO and the LUMO is less than 0.05a.u. (you will get a warning). You will also have to enter this menu for all open-shell cases other than doublets.

use file

With command use you are able to use information about occupied MOs and start vectors from a former calculation on the same molecule. file should be the path and name of the control file of this former calculation, of which all data groups related to occupation numbers and vectors will be read. As the new generated data will overwrite the existing data if both resist in the same directory, it is best and in some cases necessary to have the data of the former calculation in another directory than the one you started the define session in. Then just type use <path>/control to construct a new SCF vector from the data of the old calculation, without changing the old data. The data groups $closed shells and $open shells will be taken for your new calculation and the SCF vector from the old calculation will be projected onto the space which is spanned by your present basis set. These start vectors are usually better than the ones you could obtain by an extended Hückel calculation.

man

man allows you to declare occupation numbers or change a previous declaration manually. After selecting this command, you will get a short information about the current occupation numbers:

---------------------------------------------------------  
actual closed shell orbital selection         range  
---------------------------------------------------------  
 a1                                           #   1-  18  
 a2                                           #   1-   1  
 e                                            #   1-  13  
---------------------------------------------------------  
 
any further closed-shell orbitals to declare ? DEFAULT(y)

If you answer this question with y, you enter the orbital specification menu which will be described in Section 4.3.3.

The same procedure applies to the open-shell occupation numbers after you finished the closed-shell occupations.

hcore

hcore tells programs dscf and ridft to run without a start vector (it writes the data group $scfmo none to file control). dscf or ridft will then start from the core Hamiltonian start vector, which is the vector obtained by diagonalizing the one-electron Hamiltonian. Note that you still have to specify the occupation numbers. This concerns only the first SCF run, however, as for the following calculations the converged vector of the previous iteration will be taken. A SCF calculation with a core Hamiltonian start vector typically will take 2 - 3 iterations more than a calculation with an extended Hückel start vector (a calculation with the converged SCF vector of a previous calculation will need even less iterations, depending on how large the difference in the geometry between the two calculations is).

flip

flipping of spins at a selected atom. Requirements: converged UHF molecular orbitals and no symmetry (C1). definewill localize the orbitals, assign them to the atoms and give the user the possibility to choose atoms at which alpha-orbitals are moved to beta orbitals, or vice versa. This is useful for spin-broken start orbitals, but not for spatial symmetry breaking.

*

This command (as well as use and eht) terminates this menu, but without providing a start vector. If the keyword $scfmo exists in your input file, it will be kept unchanged (i.e. the old vector will be taken), otherwise $scfmo none will be inserted into your output file, which forces a calculation without start vector to be performed. When you leave this menu, the data groups $closed shells, $open shells (optionally) and $scfmo will be written to file. You will then reach the last of the four main menus (the General Menu) which is described in Section 4.4.

4.3.2 Assignment of Occupation Numbers

If an automatic assignment of occupation numbers is not possible or you do not except the occupation numbers generated by the EHT, you enter the following menu:

OCCUPATION NUMBER ASSIGNMENT MENU ( #e=60 #c=0 #o=0)  
 
 s         : CHOOSE UHF SINGLET OCCUPATION  
 t         : CHOOSE UHF TRIPLET OCCUPATION  
 u <int>   : CHOOSE UHF WITH <int> UNPAIRED ELECTRONS  
 l <list>  : PRINT MO’S FROM EHT IN <list>, (DEFAULT=ALL)  
 p <index> : PRINT MO-COEFFICIENTS OF SHELL <index>  
 c <list>  : CHOOSE SHELLS IN <list> TO BECOME CLOSED SHELLS  
 o <index> : CHOOSE SHELL <index> TO BECOME AN RHF OPEN SHELL  
 a <list>  : CHOOSE SHELLS IN <list> TO BECOME UHF ALPHA SHELLS  
 b <list>  : CHOOSE SHELLS IN <list> TO BECOME UHF BETA SHELLS  
 v <list>  : CHOOSE SHELLS IN <list> TO BECOME EMPTY SHELLS  
 &         : REPEAT THE EXTENDED HUECKEL CALCULATION  
 *         : SAVE OCCUPATION NUMBERS & GO TO NEXT ITEM  
 dis       : GEOMETRY DISPLAY COMMANDS  
 e         : CALCULATE EHT-ENERGY  
 f         : FURTHER ADVICE  
 <int>     = INTEGER  
 <index>   = INDEX OF MO-SHELL ACCORDING TO COMMAND s  
 <list>    = LIST OF MO-SHELL INDICES (LIKE  1-5,7-8,11)

Recommendation

Enter l to get a list of eht MO energies. Then make up your mind on what to do: closed shell, RHF open shell (not allowed for DFT) or UHF. Look at the examples below.

RHF

c 1-41,43,45 to define these levels to be doubly occupied.

UHF

a 1-5 alpha levels to be occupied, b 1-3,5 beta levels to be occupied. Or simply, s, t, or u 1 to get singlet, triplet or doublet occupation pattern.

ROHF

c 1-41,43,45 levels to be doubly occupied; o 42 level 42 should be partially occupied. You will then be asked to specify the occupation. If there are more open shells you have to repeat, since only a single open shell can be specified at a time. Watch the headline of the menu, which tells you the number of electrons assigned to MOs.

Description of Commands

s list

This command gives you a listing of all MOs and their energies as obtained from the extended Hückel calculation. For NH3 in C3v and TZVP you get, e.g.:

ORBITAL SYMMETRY  ENERGY      SHELL   CUMULATED  CL.SHL OCC.  OP.SHL OCC.  
(SHELL)   TYPE             DEGENERACY SHELL DEG. PER ORBITAL  PER ORBITAL  
    1    1a1      -15.63244      2           2      0.0000      0.0000  
    2    2a1       -0.99808      2           4      0.0000      0.0000  
    3    1e        -0.64406      4           8      0.0000      0.0000  
    4    3a1       -0.57085      2          10      0.0000      0.0000  
    5    2e         0.30375      4          14      0.0000      0.0000  
    6    4a1        0.87046      2          16      0.0000      0.0000  
TO CONTINUE, ENTER <return>

p index

This allows you to get the linear combination of basis functions which form the MO-index. Note that this refers not to the basis set you specified, but to the extended Hückel basis. index must be a single index, not an index list.

c list

This command allows you to specify closed shells. Their occupation will be 2 per MO, the total occupation the shell degeneracy which you can obtain by using command s. list is a list of shell indices like 1-13 or 1,3-5,7.

o index

This command allows you to specify open shells. index must be a single shell index, not an index list. You will then be asked for the number of electrons per MO which shall be contained in this shell. For example, for a fluorine atom you should choose o n (where n is the index of the p-shell) and an occupation of 53 per MO. You may enter the occupation numbers as simple integers or as integer fractions, e.g. 1 for the s-occupation in sodium, 5/3 for the p-occupation in fluorine.

v list

With this command you can remove an orbital occupation, if you specified a wrong one. list is again a list of shell indices in usual syntax.

&

This command has a different meaning in this menu than in the rest of define. Here it will repeat the extended Hückel calculation (perhaps you want to change some Hückel parameters for the next one).

*

* will not bring you back to the occupation numbers menu, but will terminate the whole occupation number and start vector section and will bring you to the last main menu, which is described in Section 4.4. If you want to leave this menu without assigning all electrons in your molecule to shells, define will issue a warning and suggest to continue defining occupation numbers. You can ignore this warning, if you do not want to assign all electrons.

e

Calculates and displays the extended Hückel total energy of your molecule.

f

f will give you some information about the commands in this menu.

You may overwrite occupation numbers once given by just redefining the corresponding shell. For example, if you choose shells 1–10 as closed shells and afterwards shell no. 9 as open shell (with any occupation number), the open shell will be correctly assigned.

4.3.3 Orbital Specification Menu

define provides the possibility to assign the occupation numbers of the MOs manually, if you like. To do that, use the command man in the occupation number main menu and you will arrive at the following submenu:

------------- ORBITAL SPECIFICATION MENU --------------  
 
<label> <list>  : select orbitals within <list>  
-<label> <list> : skip orbitals within <list>  
&               : ignore input for last label  
clear           : clear all assignments  
p(rint)         : print actual orbital selection  
for help, type ? or help // for quit, type * or q(uit)


Depending on whether you are in the closed- or in the open-shell section, the commands of this menu refer only to the corresponding type of orbitals. The commands of this menu do not need much explanation. <label> is the irrep label of one irreducible representation of the molecular point group (e.g. a1, b2, t1g, …). <list> is a list of orbital indices within this irrep (e.g. 1,2,4 or 1-8,10,11). p or print will give you the same listing of the orbital occupations as you saw before entering this menu. After you leave this submenu, you will be back in the occupation numbers main menu.

4.3.4 Roothaan Parameters

In open-shell calculations within the restricted Hartree–Fock ansatz (ROHF), the coupling between the closed and the open shells must be specified using two parameters a and b, which depend on the type of the open shell, the number of electrons in it (the electron configuration), but also on the state to be calculated. For example, there are three states arising from the s2p2 configuration of an atom (3P, 1D, 1S) which have different values of a and b. For the definition of these parameters and their use refer to Roothaan’s original paper [26]. For simple cases, define sets these parameters automatically. If not, you have to enter them yourself. In this case, you will get the following message:

ROOTHAAN PARAMETERS a AND b COULD NOT BE PROVIDED ...  
TYPE IN ROOTHAAN a AND b AS INTEGER FRACTIONS  
OR ENTER val FOR AN AVERAGE OF STATES CALCULATION  
OR ENTER  &  TO REPEAT OCCUPATION NUMBER ASSIGNMENT

Note that not all open shell systems can be handled in this way. It is possible to specify a and b for atomic calculations with sn,pn,d1, and d9 configurations and for calculations on linear molecules with πn and δn configurations. Furthermore, it is possible to do calculations on systems with half-filled shells (where a=1, b=2). In the literature you may find tabulated values for individual states arising from dn configurations, but these are not correct. Instead, these are parameters for an average of all states arising from these configurations. You can obtain these values if you enter val on the above question. For a detailed description see Section 6.3.

4.3.5 Start-MOs for broken symmetry treatments ("flip")

Broken-symmetry treatments suggested by e.g. Noodleman or Ruiz are a popular tool for the calculation of spin coupling parameters in the framework of DFT. As an example one might consider two coupled CuII centers, e.g. for a (hypothetical) arrangement like this:

$coord  
 0.0  2.7 0.0 cu  
 0.0 -2.7 0.0 cu  
 0.0 -6.1 0.0 f  
 0.0  6.1 0.0 f  
 2.4  0.0 0.0 f  
-2.4  0.0 0.0 f  
$end

The high-spin case, a doublet with an excess alpha electron at each Cu atom, "aa" in an obvious notation, preserves D2h symmetry, while the low spin state "ba" does not. For a broken-symmetry treatment it is advidsable to calculate the high-spin state first, and then broken-symmetry state(s); from the energy difference(s) one may calculate approximate values for the spin-spin coupling parameters as described by e.g. the above authors. Access to broken-symmetry states usually is possible by the choice of appropriate start-MOs, followed by an SCF-procedure. Start MOs may be obtained by first applying a localization procedure to the MOs of the high-spin state and then by "moving" localized alpha orbitals to the beta subset.

The preparation of broken-symmetry start-MOs can be done with define (semi-)automatically. Prerequisite is a converged wave function for the high-spin state in C1-symmetry, that fulfills the aufbau principle.

If in this case one enters flip in the orbital definition menu, define selects the occupied valence orbitals of the system (by an orbital energy criterion, which one can usually accept, unless the system is highly charged and the orbital energies are shifted). Next a Boys orbital localization procedure is carried out, which - depending on the size of the problem - may take some time. Then the user is asked:

ENTER INDICES OF ATOMS OR ELEMENT TO BE MANIPULATED (example: 1,2-3 or "mn")

In case of our above example one may enter "cu", which immediately leads to the following output (a def-SV(P) basis and the B-P functional were used for the high-spin state):

 RELEVANT LMOS FOR ATOM   1 cu  
 ALPHA:  
index occupation "energy"     s    p    d    f   (dxx  dyy  dzz  dxy  dxz  dyz)  
  15    1.000     -0.357     0.01 0.00 0.98       0.20 0.27 0.01 0.50 0.00 0.00  
  18    1.000     -0.357     0.01 0.00 0.98       0.20 0.27 0.01 0.50 0.00 0.00  
  20    1.000     -0.335     0.00 0.00 1.00       0.00 0.00 0.00 0.00 1.00 0.00  
  22    1.000     -0.333     0.01 0.00 0.99       0.13 0.03 0.32 0.00 0.00 0.51  
  23    1.000     -0.333     0.01 0.00 0.99       0.14 0.03 0.34 0.00 0.00 0.49  
 
 BETA:  
  39    1.000     -0.326     0.00 0.00 1.00       0.33 0.08 0.09 0.00 0.00 0.50  
  41    1.000     -0.326     0.00 0.00 1.00       0.33 0.08 0.09 0.00 0.00 0.50  
  43    1.000     -0.321     0.00 0.00 1.00       0.00 0.00 0.00 0.00 1.00 0.00  
  46    1.001     -0.318     0.05 0.00 0.95       0.00 0.43 0.51 0.00 0.00 0.00  
 
 RELEVANT LMOS FOR ATOM   2 cu  
 ALPHA:  
index occupation "energy"     s    p    d    f   (dxx  dyy  dzz  dxy  dxz  dyz)  
  16    1.000     -0.357     0.01 0.00 0.98       0.20 0.27 0.01 0.50 0.00 0.00  
  17    1.000     -0.357     0.01 0.00 0.98       0.20 0.27 0.01 0.50 0.00 0.00  
  19    1.000     -0.335     0.00 0.00 1.00       0.00 0.00 0.00 0.00 1.00 0.00  
  21    1.000     -0.333     0.01 0.00 0.99       0.13 0.03 0.32 0.00 0.00 0.51  
  24    1.000     -0.333     0.01 0.00 0.99       0.14 0.03 0.34 0.00 0.00 0.49  
 
 BETA:  
  40    1.000     -0.326     0.00 0.00 1.00       0.33 0.08 0.09 0.00 0.00 0.50  
  42    1.000     -0.326     0.00 0.00 1.00       0.33 0.08 0.09 0.00 0.00 0.50  
  44    1.000     -0.321     0.00 0.00 1.00       0.00 0.00 0.00 0.00 1.00 0.00  
  45    1.001     -0.318     0.05 0.00 0.95       0.00 0.43 0.51 0.00 0.00 0.00  
 
 a2b    : FLIPPING ALPHA TO BETA (default)  
 b2a    : FLIPPING BETA TO ALPHA  
 r      : repeat atom choice  

As evident from the second column, for each Cu atom five localized alpha and four localized beta orbitals were generated which are of d-type (the sixth column labelled "d" shows values close to 1, the other columns such close to 0). The six columns at the right show the individual contributions of the six cartesian d-functions.

What has to be done to generate start MOs for the "ba"-case? Obviously one of the five localized alpha spin orbitals from the first Cu atom (atom label 1 cu) has to become a beta spin orbital. These five orbitals have the indices 15, 18, 20, 22, 23. In order to avoid linear dependencies, it is advisable to take the orbital that has no beta-analogue. This can be found by comparing the contributions of the six d-functions. In the present example this is the case for the localized alpha orbitals 15 and 18: in contrast to all localized beta orbitals they show significant contributions from dxy. One thus enters

a2b  
15

and after confirming the replacement of original MOs with the generated start-MOs one is finally asked

 It is advisable to modify damping and orbital shift in the following way:  
 
 $scfdamp   start=5.000  step=0.050  min=0.500  
 $scforbitalshift  automatic=1.0  
 $scfiterlimit 999  
 
 Do you want to replace the corresponding entries in the control-file? (y)

which should be confirmed, as otherwise the prepared spin state might be destroyed during the SCF iterations.

From this input one may start the SCF(HF/DFT)-procedure. For recommended choices of DFT functionals and formulae to calculate the coupling parameters from these energy differences please consult the papers of the above-mentioned authors. For reasons of economy, a pre-optimization by a pure (non-hybrid) DFT-functional is reasonable.

Important: For the converged wave function one should check, whether the resulting state is really the desired one. This can quite reliably be done by a Mulliken population analysis. For this purpose, add $pop to the control file, type ridft -proper or dscf -proper, respectively, and check the signs of the calculated numbers of unpaired electrons in the output.