The Atlas ATMOL Manual is no longer available. This version of the Manual dates from 1982 after ATMOL was ported to the Cyber 205 at UMRCC.
In the following text, a program capable of performing closed and open shell restricted and unrestricted Hartree-Fock calculations and Boys localisation of molecular orbitals, is described.
The modules available are as follows:
In default the program will request 3 large pages of main memory, this can be varied by means of the pre-directives LPAGE or MEMORY [1]. All pre-directives [1] are applicable and should be presented before the program specific directives described below.
Data input and printed output are on FORTRAN streams 5 and 6 respectively. Punched output is sent to FORTRAN stream 7. These streams need not be explicitly mentioned in the JCL.
The following data sets will be used by the SCF program, and should be mentioned in the JCL, in REQUEST, ATTACH, MFLINK or GETFEP commands:
MAIN FILE: this contains the atomic 2-electron integrals generated via the integrals program [2]. This data set can be assigned any AFN file ED0 to ED6 or MT0 to MT7, the default is ED2.
DUMP FILE: this data set contains general information such as the geometry and basis set of the system under investigation as well as the atomic 1-electron integrals. The SCF program also outputs the eigen vectors on termination of the SCF program to the DUMP FILE. This data set can be assigned any AFN file ED0 to ED6 or MT0 to MT7, the default is ED3.
SCRATCH FILE: the SCF program requires a data set to be made available to be used as a scratch area for transient data handling. This data set is assigned to the AFN file ED7.
Data is read from FORTRAN stream 5. Th syntax of the first line is:
PROG NBASIS NCORE IBLKD [ DDUMP [ INTSEC ] ]
where
The remainder of the data consists of directives. The directive structure applicable to each program will be outlined in this manual.
The directives for the closed shell SCF program may be presented in any order, unless otherwise specified.
This directive allows the user to define an 80 character title for the run, and extends over two data lines. The first line contains the character string TITLE in the first data field, the second line contains the required 80 character title.
Example :
TITLE
H2O CLOSED SHELL
This directive consists of one data line, read to variables TEXT, MAXC using format (A,I).
TEXT should be set to the character string MAXCYC.
MAXC is an integer used to specify the maximum number of iteration cycles required.
The directive may be omitted when MAXC will be set to the default value of 50.
The following conditions cause termination of iteration:
This directive is used to define the location of the MAIN FILE to be used as a source of 2-electron integrals, and consists of four data lines, the first containing the character string MFILE in the first data field.
Consider the case where the MAIN FILE is split into LTAPE sections, then the remaining 3 lines of the MFILE directive is consists of LTAPE data fields.
The second line is read to variables (DDMAIN(I), I=1,LTAPE) in A-format. DDMAIN(I) specifies the ATMOL file used to assign the data set in which the I'th section of the MAIN FILE is to be found. The value of LTAPE is determined by the program from the number of data fields on this line.
The third line is read to variables (IBLK(I), I=1,LTAPE) in I-format. IBLK(I) specifies the starting block number of the I'th section of the MAIN FILE.
The fouth line is read to variables (LBLK(I), I=1,LTAPE) in I-format. LBLK(I) defines the block number of the first vacant or unused block after the I'th section of the MAIN FILE. Thus the number of blocks in the I'th section will be given by LBLK(I)-IBLK(I). If the I'th section is known to be an ENDFILE block, LBLK(I) may be set to zero.
Example :
A MAIN FILE consists of three sections, stored on data sets assigned as ED2,MT2,MT3 all starting at block 1, and all ending with ENDFILE blocks. The appropriate MFILE directive is as follows:
MFILE
ED2 MT2 MT3
1 1 1
0 0 0
This directive consists of a single data line, read to variables TEXT, ISIZE using format (A,I).
TEXT should be set to the character string SIZE.
ISIZE is an integer used to specify the maximum size (in blocks) to which the DUMP FILE is allowed to grow.
The SIZE directive may be used to overide the previous SIZE assignment as given by the integrals program.
Example :
SIZE 90
The CTRANS directive may be used to define a basis set for the SCF calculation, where the basis functions comprise linear combinations of the orginal basis as defined at the integral evaluation [2]. The acronym LCBF (Linear Combined Basis Function) will be used when referencing a linear combination defined by the CTRANS directive. The first data line is read to variables TEXT,N using format (A,I).
TEXT should be set to the character string CTRANS.
N is an integer specifying that the first N 'LCBF' are to equal to the first N 'original' basis functions. N may be omitted, when it is given value zero.
Data concerning further LBCF to be appended to the list generated by the first directive is presented next. Each LCBF is introduced in turn. If a LCBF required NTRAN terms, NTRAN data lines would be required for its definition. The first being read to variables C,ITRAN,NTRAN using format (F,2I). The ITRAN'th 'original' basis function will be included in the LCBF with coefficient C. The value of NTRAN specifies the number of terms within the linear combination, and NTRAN-1 successive data lines will then be read, (if NTRAN=1, no further lines required), each being read to variables C,ITRAN using format (F,I). The ITRAN'th 'original' basis function will be included with coefficient C. When all data concerning all the LCBF has been supplied, the directive is terminated by a line containing the character string END in the first data field.
The following points may prove helpful:
Example 1:
An integral evaluation was performed with 8 basis functions. The calculation is to be performed in the following basis set: LCBF 1 to 6 are to be the same as 'original' basis functions 1 to 6, and symmetric and anti-symmetric combinations of the 'original' basis functions 7 and 8 are to be taken. The CTRANS data would be presented as:
CTRANS 6
1 7 2
1 8
1 7 2
-1 8
END
Example 2:
An integral evaluation was performed with 15 basis functions. It is necessary to perform a calculation with 'original' basis functions 11 to 15 removed. The CTRANS data would consist of:
CTRANS 10
END
For small cases, the use of the CTRANS directive to 'delete' basis functions is convenient. For larger basis sets, the MERGE facilities of the INTEGRAL [2] and SERVICE [3] programs will be found to save SBU time.
Example 3:
An integral evaluation run has been performed using S,P and D groups of basis functions, the COMBINE option (INTEGV module only) not being used used. The 'original' basis set is therefore, in the order, s,p(x),p(y), p(z),d(xy),d(xz),d(yz),d(x**2),d(y**2) and d(z**2). If the user requires linear combinations of the D basis functions so that the linear combinations span only the D irreducible representation of the spherical point group, as compared with the 'original' D basis, which spans the S and D irreducible representations. The CTRANS directive may be used as follows:
CTRANS 7
1 8 2
-1 9
-1 8 3
-1 9
2 10
END
This directive is used to control the quantity of printed output produced at convergence (or at non-convergence if the OUTPUT directive is specified). The FPRINT directive may be used to increase the volume of printing, which consists of one data line set to the character string FPRINT in the first data field. Subsequent data fields are read in A-format, and specify that printing as detailed in Table 1 is to occur at convergence.
Table 1: Paramaters of the FPRINT Directive
___________________________________________
Data Field Printing of
__________ ___________
S Matrix of overlap integrals over the
LCBF.
T Matrix of kinetic energy integrals
over the LCBF.
T+V Matrix of integrals for the 1-electron
part of the Hamiltonian operator.
FOCK Roothaan-Hartree-Fock matrix.
POPS Basis function Mulliken overlap
populations.
R The 1-particle density matrix.
X The matrix of integrals over the X
component of the dipole operator.
Y The matrix of integrals over the Y
component of the dipole operator.
Z The matrix of integrals over the Z
component of the dipole operator.
Example :
FPRINT FOCK POPS
will cause printing of the Fock and overlap population matrices at convergence (or non-convergence if the OUTPUT directive is included).
This directive consists of one data line, the first data field of which contains the character string PUNCH. Subsequent data fields are in A-format, and specify that 'punching' (to FORTRAN stream 7) is to occur at convergence (or non-convergence if the OUTPUT directive with a CARDS paramater is used). The other data fields may be set to:
VECTORS specifies that converged eigen vectors are to be punched.
MOPOPS specifies that the population numbers of the converged molecular orbitals are to be punched.
EVALUES specifies that the converged eigen values are to be punched.
If the second paramter is omitted and the character string PUNCH is included, this is equivalent to the coding
PUNCH VECTORS
The format used for the punched output:
VECTORS: The first line punched will be the title of the run (as read from the TITLE directive). The second punched line contains the text VECTORS in columns 1 to 7, a space character in column 8, and in column 9 to 16 the FORTRAN format to be used for the punching of the vectors, as specified by the FORM1 parameter of the FORMAT directive, or to 5F16.9 in the absence of the FORMAT directive. The eigen vectors are then punched by means of code of the form:
WRITE(7,FORM1)((Q(I,J),I=1,NEWBAS),J=1,NEWBAS)
where FORM1 denotes the format specified by the FORM1 parameter of the FORMAT directive (or default 5F16.9).
MOPOPS: The first line punched consists of the title of the run. The second line contains the text MOPOPS in columns 1 to 6, a space character in column 7, and the FORTRAN format (as specified by the FORM2 parameter of the FORMAT directive, or 5F16.9 in default) to be used for punching the molecular orbital populations is punched in columns 8 to 15. The molecular orbital populations are then punched, by means of the code:
WRITE(7,FORM2)(POP(I),I=1,NEWBAS)
In the case of the closed shell SCF program, the first NCORE populations will equal 2, the remainder being zero.
EVALUES: The first line punched consists of the title of the run. The second line contains the text EVALUES in columns 1 to 7, followed by a space character in column 8, an then the FORTRAN format (as specified by the FORM3 parameter of the FORMAT directive or 5F16.9 in default) to be used for punching the eigen values in columns 9 to 16. The eigen value are then punched by means of the code:
WRITE(7,FORM3)(EIG(I),I=1,NEWBAS)
This directive is used for specifying the FORTRAN formats to be used for punching VECTORS,MOPOPS and EVALUES. The directive consists of one data line read to variables TEXT,FORM1,FORM2,FORM3 using format (4A).
TEXT should be set to the character string FORMAT.
FORM1,FORM2,FORM3 are the parameters referred to in the description of the PUNCH directive which control the printing of the VECTORS, MOPOPS and EVALUES.
If the paramter of the FORMAT directive is omitted, it will be given the default value 5F16.9.
Example :
FORMAT 4D20.11 4F20.8 4D20.9
will cause the punching of VECTORS,MOPOPS and EVALUES to be in FORTRAN formats 4D20.11,4F20.8 and 4D20.9 respectively.
This directive consists of one data line, read to variables TEXT, TEST using format (2A).
TEXT should be set to the character string OUTPUT.
TEST may optionally be set to the character string CARDS, or may be omitted.
The OUTPUT directive is only of relevance if non-convergence occurs, and will then initiate a full printout of the current state of the calculation. If the parameter CARDS appears on the OUTPUT directive, any punching requested by the PUNCH directive will also be produced.
Omission of the OUTPUT directive causes the program to terminate execution without printing or punching in a non-convergence situation.
This directive consists of a single data line with the character string IPRINT in the first data field. If used causes printing of the orbital energies resulting from each iterative cycle.
This directive consists of one data line, read to variables TEXT,B,K using format (A,F,I).
TEXT should be set to the character string ACCURACY.
B,K a convergence threshold, T = B*10**(-K) will be computed. At convergence, the elements of the density matrix will be converged to within an absolute error T.
The ACCURACY directive may be omitted, then the default value T=10**(-5) will be used.
Example 1:
ACCURACY 1 4
sets the convergence threshold to 10**(-4).
Example 2:
ACCURACY 1 9
will cause the convergence threshold to be set to 10**(-9), this being its lowest possible value.
Example 3:
ACCURACY 5 7
will cause the convergence threshold to be set to 5*10**(-7).
This directive may be used to set the maximum SBU time allowed for the calculation, and consists of one data line read to variables TEXT, TIMEMX using format (A,F).
TEXT should be set to the character string TIME.
TIMEMX is used to specify the time maximum in units of SBUs. The program will terminate iterative cycling when the SBU time use for the job, plus the iteration cycle time is greater than TIMEMX.
The TIME directive is normally omitted, when TIMEMX will be set by the program, to the number of SBU's of total job time as requested on the user's job card or the time remaining after a sequence of job steps.
Example :
TIME 750
causes the program to terminate cycling, whether convergence has been achieved or not, when it's found that
(total jobtime used + cycle time) > 750 SBUs.
Because of the 'level shifting' technique for forcing convergence, the possibility arises that convergence to a stationary point on the energy surface may occur, such that the converged molecular orbitals do not represent the ground state. To circumvent this problem, the program will automatically, if required, switch virtual and occupied orbitals so as to force convergence to a configuration which obeys the 'aufbau' principle. Such switching is normally allowed to occur only when the presumed occupancy of the input molecular orbitals is within a tolerance of 0.075.
The AUTO directive may be used to override these usual conventions, and consists of a single data line, read to varaibles TEXT, SWITCH using format (A,F).
TEXT should be set to the character string AUTO.
SWITCH automatic orbital switching will not be allowed until the molecular orbitals of a non-aufbau configuration achieve an accuracy within a tolerance specified by SWITCH.
The following notes may prove helpful to the user employing this directive:
Example 1:
The automatic orbital switching may be completely suppressed by the directive of the form:
AUTO 0
This may be useful in the study of excited closed shell configurations which do not obey the aufbau principle.
Example 2:
An optimal convergence rate may often be achieved by allowing automatic orbital switching from an early stage in the calculation.
AUTO 100
Will permit switching to occur from the start of the calculation. The user is asked to note that the strength of the convergence guarantee provided by the level shifting feature is weakened if switching is used, so that convergence problems may arise if switching is allowed at a too early stage in the calculation. Such a situation may arise when the default value of SWITCH is invoked, although this is not likely, and may be remedied by the reduction in the value of SWITCH, by means of a directive:
AUTO 0.001
This consists of one data line read to variables TEXT,D1,E1,IBRK,D2,E2 using format (A,2F,I,2F).
TEXT should be set to the character string DAMP.
D1 is the damp factor up to iterative cycle specified by IBRK.
E1 is the level shifter up to iterative cycle specified by IBRK.
IBRK is an integer used to a specify the cycle number.
D2 is the damp factor after the iterative cycle given by IBRK.
E1 is the level shifter after the iterative cycle given by IBRK.
An alternative form of DAMP is permitted, consisting of only three parameters, read to variables TEXT,D1,E1 using format (A,2F). If used, this set the IBRK parameter to the default 999, D2 and E2 will be given the default values of 1 and 0 respectively. D1 and E1 have there usual meanings.
The DAMP directive may be omitted, then the program will assign the following default settings:
D1,D2=1.0 E1,E2=0.0 and IBRK=999
The choice of the above defaults corresponds to performing the iterations in precisely the manner prescribed by Roothaan [4].
A full treatment of the role of damp factors and level shifters will be discussed later. At present one shall confine the discussion to the practical matter of selecting reasonable parameters. Experience indicates taht, at least for ground state studies, there is little to be gained by choosing values for the damp factor (D1,D2) different from unity, and further discussion will assume that the value of unity has been chosen. It remains to consider the values of (E1,E2), level shift parameters.
Let K denote a parameter which notionally measures the rate of convergence of the iterations. When K is positive, the procedure is convergent, the larger the value of K, the more rapid the convergence. When K is negative, the process is divergent.
A common feature of all convergence curves is that convergence sets in at some critical value of the level shifter, the optimum rate of convergence being found at somewhat higher values of level shifter. Use of a level shifter with a value greater than the optimum does not prevent convergence, but convergence is slowed down, whilst use of a level shifter which is too small may cause divergence.
The user may be helped in the choice of the level shifter parameters if they take notice of the following points:
DAMP 1 1 3 1 0
will often prove successful.For further reading see also Damp Factors and Level Shifters for the Closed Shell SCF Program.
Seven directives are available for obtaining trial eigen vectors to act as input to the first iterative cycle, namely VECTORS,START,ALPHAS,GETQ, RESTORE,EXTRA and COMBINE. The generic term LOADQ is not itself a directive. One, and only one, LOADQ directive can appear in the input data, the structure of these directives are now discussed.
This directive consists of one data line read to variables TEXT, FORMV using format (2A).
TEXT should be set to the character string VECTORS.
FORMV is used to specify the FORTRAN format to be used to input trial eigen vectors, or may be omitted, when the default is 5F16.9.
The remainder of the data for the VECTORS directive is read using the code:
READ(5,FORMV)((Q(I,J),I=1,NEWBAS),J=1,NEWBAS)
where FORMV denotes the format specified by the FORMV parameter to be found on the directive initiator. Vectors punched out by the SCF programs will be found suitable for re-input, so long as the user remembers to remove the first card punched, which contains the TITLE of the job which was responsible for producing the vectors.
This consists of one data line with the character string START in the first data field. Trial vectors will be formed by diagonalisation of the 'unscreened' 1-electron hamiltonian operator matrix. The vectors will be ordered according to increasing eigen value, and the NCORE vectors of lowest eigen value deemed to constitute the trial set of occupied orbitals.
The first data line consists of the character string ALPHAS in the first data field. Subsequent lines should contain NEWBAS real numbers read in free F-format to a vector:
(A(I),I=1,NEWBAS)
Hence the sequences:
ALPHAS
5.5 2.1 .5 .5 .3
and
ALPHAS
5.5
2.1
.5
.5
.3
are equivalent. A(I) should be set equal to the negative of the expected value of the I'th element of the Roothaan- Hartree- Fock martrix (in basis set representation) at self consistency. The utility of the ALPHAS directive stems from the fact that diagonal Roothaan- Hartree-Fock matrix elements remain approximately invariant under change in molecular environment for a given basis set. For example, a study of the water yields ALPHAS parameters appropriate for the oxygen and hydrogen atoms in other molecular environments in the basis set used. The user, if applying the same basis set, can quickly build a 'library' of appropriate ALPHAS parameters, simply by printing out converged Roothaan-Hartree-Fock matrices by means of the FPRINT directive:
FPRINT FOCK
If a reasonably accurate set of ALPHAS parameters is available, the ALPHAS directive will provide a good mechanism for the trial eigen vectors. If the user is unable to provide good estimates for the ALPHAS parameters, the START directive may be used.
The following algorithm is used in the implementation of the ALPHAS directive. The matrix representatives of the kinetic energy operator (T) and nuclear attraction operator (V) are internally available to the program. An approximation for the coulomb/exchange 2-electron operator (G) and an approximate Fock operator (H) could be constructed from:
H = T + V + G
The user supplies the A-parameters, which are approximations to the negative to the diagonal Fock operator elements.
Hii = -A(I)
The diagonal elements of G can therefore be formed from:
Gii = Hii - Tii - Vii
where the guessed values of Hii are used.
Assume that the off-diagonal elements of G and V can be expressed in the form of Mulliken approximations:
Gij = Kij * (Gii + Gjj) Vij = Kij * (Vii + Gjj)
where Kij is some unknown scalar. Elimination of Kij gives the result:
Gij = Vij * (Gii + Gjj) / (Vii + Vjj)
whence approximations to the off-daigonal elements of the Fock operator can be computed from:
Hij = Tij + Vij + Vij * (Gii + Gjj) / (Vii + Vjj)
The eigen value problem:
HQ = SQA
is solved with Fock operator (H) approximated as above. The resulting eigen values (Q) are arranged in order of ascending eigen values, and define the trial molecular orbitals as linear combinations of the LCBF.
This directive consists of a single data line, read to variables TEXT, ISECT using format (A,I).
TEXT should be set to the character string RESTORE.
ISECT is an integer used to specify the section of the DUMP FILE where trial eigen vectors may be found. Such trial vectors will have normally been placed on the DUMP FILE by a previous run of the SCF program, the section number used to store the vectors being specified by the ENTER directive of a previous run.
Example :
RESTORE 1
Eigenvectors are restored from section 1 of the DUMP FILE, such vectors having been routed to section 1 by means of a previous run of the SCF program where the directive
ENTER 1
was used.
This directive causes restoration of eigen vectors from a 'foreign' DUMP FILE, and consists of a single data line read to variables TEXT, DDFORN, IBLKF, ISECTF using format (2A,2I).
TEXT should be set to the character string GETQ.
DDFORN is the ATMOL file name assigned to the 'foreign' DUMP FILE.
IBLKF is an integer specifying the starting block number of the 'foreign' DUMP FILE.
ISECTF is an integer used to specify the section wherein the required vectors are to be found on the 'foreign' DUMP FILE.
Example :
GETQ ED5 1 3
This example indicates the recovery of the eigen vectors from section 3 of the data set assigned to ED5 starting at block 1.
If the present set of LCBF's differs from that used in a previous calculation by the addition of a number of LCBF's. Also that these additional LCBF's do not constitute a significant influence in the newer calculation. Then molecular orbitals of the previous calculation would provide a good starting point for the new calculation, thus the EXTRA directive has been implemented to allow the previous vectors to be used as starting vectors. The first data line is read to variables TEXT, DDEXT, IBLKE, ISECTE using format (2A,2I).
TEXT should be set to the character string EXTRA.
DDEXT specifies the ATMOL file name assigned to the data set, where the trial vectors reside (in the smaller basis set).
IBLKE is an integer specifying the starting block number of the foreign data set.
ISECTE is an integer specifying the section number where the trial vectors are to be found. This integer was specified on the ENTER directive of the run which created the trial vectors.
If the number of LCBF in the reference calculation is denoted by OLDBAS. Then the number of additional LCBF required in the new calculation is given by NEWBAS-OLDBAS=NEXTRA. Subsequent data lines of the EXTRA directive will consist of NEXTRA integers read in free I-format, which are used to specify the LCBF which appear in the new calculation, and which do not appear in the reference calculation.
The last data field should contain the character string END. This entry will issue the ORTHOG directive in SCHMIT mode so that the trial eigen vectors will be automatically Schmit orthogonalised.
Example 1:
EXTRA ED5 1 1
11 12 13 14 15 END
This example specifies that the trial vectors reside on ED5 starting at block 1 of section 1. LCBF 11 to 15 appear in the new calculation but do not appear in the reference calculation. It is assumed that those basis functions of the new calculation in common with the reference calculation occur in the same order. The above example may be specified in several ways:
EXTRA ED5 1 1
11
12
13
14
15
END
and
EXTRA ED5 1 1
11 TO 15
END
are all equivalent. In the last permutation the character string TO was employed to link a sequence of consecutively numbered additional LCBF.
Example 2:
The sequences:
EXTRA ED6 200 2
1 5 6 7 8 20 21 22 23 END
and
EXTRA ED6 200 2
1 5 TO 8 20 TO 23 END
are equivalent.
If the user has performed separate SCF calcuations on two systems, A and B, the converged vectors being stored on their DUMP FILEs. And if the user now wishes to perform a calculation on the system AB, where A is sited at a considerable distance from B. The trial vectors from A and B, if combined would provide a suitable starting point for the calculation on AB. The COMBINE directive has been supplied to provide such a facility. The first data line is read to variables TEXT, DDA, IBLKA, ISECTA, DDB, IBLKB, ISECTB using format (2A,2I,A,2I).
TEXT should be set to the character string COMBINE.
DDA,IBLKA,ISECTA the 'A' set of vectors will be read from section ISECTA of a DUMP FILE commencing at block IBLKA on an ATMOL file specified by the parameter DDA.
DDB,IBLKB,ISECTB the 'B' set of vectors will be read from section ISECTB of a DUMP FILE commencing at block IBLKB on an ATMOL file specified by the parameter DDB.
In subsequent data lines the 'origins' of the LCBF of the system AB are defined. Such data lines read to variables I,J,AORB,K using format (2I,A,I). This means that LCBF I to J (inclusive) of the AB calculation are in one to one correspondence with the LCBF K onwards of either the A or B calculation, depending upon the specification of AORB. AORB should be set to the character A or B accordingly. When all the AB LCBF have had their 'origins' specified, this part of the data input is terminated by the character string END in the first data field.
Note that all the LCBF of the system AB must be in correspondence with LCBF of either A or B, and that the total number of LCBF in the A and B systems must equal the number in the AB system.
The next part of the COMBINE data specifies how the trial eigen vectors for the AB system are derived. These data lines are read to variables I,J,AORB,K using format (2I,A,I), and mean that trial vectors I to J (inclusive) of the AB system are to be in one to one correspondence with vectors K onwards of either the A or B set of vectors, according to the specification of AORB. AORB should be set to the character A or B. When all the AB trial vectors have been defined, the directive is terminated with a data line containing the character string END in the first data field. The directive will issue an ORTHOG directive in SCHMIDT mode.
Example :
LCBF 1 to 5 of AB are to be in one to one correspondence with LCBF 1 to 5 of A. LCBF 6 to 10 of AB correspond to LCBF 1 to 5 of B. LCBF 11 to 20 of AB corresponds to LCBF 6 to 15 of A. LCBF 21 to 26 of AB are to be in correspondence with LCBF 6 to 11 of B. The trial vector matrix of AB is constructed as follows:
Vectors 1 to 4 are to correspond to vectors 1 to 4 of the A system. Vectors 5 to 7 are to correspond to vectors 1 to 3 of the B system. Vectors 8 to 18 are to correspond to vectors 5 to 15 of the A system. Vectors 19 to 26 are to corespond to vectors 4 to 11 of the B system.
The COMBINE data would look like:
COMBINE ED4 1 1 ED4 200 3
1 5 A 1
6 10 B 1
11 20 A 6
21 26 B 6
END
1 4 A 1
5 7 B 1
8 18 A 5
19 26 B 4
END
where the 'A' set of vectors are retrieved from section 1 starting at block 1 on a DUMP FILE assigned to ED4. While the 'B' set of vectors are stored on section 3 starting at block 200 of the DUMP FILE assigned to ED4.
The following notes on the LOADQ directives may prove helpful to the user:
This directive consists of a single data line set to the character string SDIAG in the first data field. This directive issues a command to output the overlap between the basis functions.
It may be used to indicate if by accident in the integral generation, a basis function which has been declared more then once on the same centre, this causes the SCF program to abort.
This directive is used to cause orthonormalisation of the trial eigen vectors, and consists of one data line read to variables TEXT,TEST using format (2A).
TEXT should be set to the character string ORTHOG.
TEST may be set to one of the character strings SCHMIDT or SYMM or may be omitted. If SYMM is specified, the Lowdin symmetric S**(-0.5) orthormalisation scheme will be used. If TEST is omitted, or if SCHMIDT is specified, then the Gramm-Schmidt orthonormalisation scheme will be used.
Note if the directive is issued twice, once with the SCHMIDT parameter, and once with the SYMM parameter, then the symmetric scheme will be chosen. Note also the EXTRA and COMBINE directives issue ORTHOG directives internally.
This directive is used to cause no-orthonormalisation of the trial eigen vectors, and consists of one data line set to the character string NOORTH.
This directive turns off the SCHMIDT and SYMM orthonormalisation flags. The use of this directive may be seen if no-orthonormalistion is wished after the use of the LOADQ directives of EXTRA and COMBINE. This directive must be used after either one of these LOADQ directives, if to invoke no-orthonormalisation.
The first data line should contain the character string SWAP in the first data field. Subsequent data lines are read to variables I,J using format (2I).
The effect is that the I'th and J'th molecular orbitals, as entered into by the LOADQ directive, are interchanged. When all the interchanging lines have been presented, the directive should be terminated by a data line containing the character string END in the first data field. The following notes may prove helpful:
Example :
SWAP
10 12
11 13
END
Molecular orbitals 10 and 11 are interchanged with molecular orbitals 12 and 13 respectively.
In order to make effective use of this directive, the basis set should be symmetry adapted (CTRANS directive). The purpose of the directive, which must appear after the LOADQ directive, is to allow reordering of the trial molecular orbitals according to their symmetry labels. The first data line contains the character string ORBSYM in the first data field.
Following this text should be the basis function symmetry declaration lines, at least one line per each irreducible representation is required. The first data field of such a line is read in A-format, and should contain a user specified symmetry species label. Subsequent data fields are read in I-format. Let the value of the integer specified in such an I-field be j. Then the j'th basis function will be assigned the symmetry label found in the first data field of the line. The following:
SIGMAG 1 2 3 4 7 8 9 10
comprises a valid basis function symmetry designation line. Such a line may be shortened if sequences of consecutive integers are presented, by means of the character string TO. Thus the above list may be expressed:
SIGMAG 1 TO 4 7 TO 10
The basis function symmetry designation lines are presented until all basis functions have been given symmetry labels, and a data line which contains the character string END in the first data field is used to terminate the sequence. Notice that it is therefore impossible to use the text END as a symmetry label. Any basis functions not assigned a symmetry label by the user will be assigned the label * by the program. From the symmetry labels of the basis functions, the program works out the symmetry labels of the trial molecular orbitals.
The next phase of the ORBSYM directive consists of one data line whose purpose is to specify the ordering of the trial molecular orbitals according to symmetry label. The first data field is read in A-format, and should be equal to a symmetry label, as defined by one of the basis function symmetry designation lines presented above (or possibly '*'). The remaining data fields are read in I-format. Let the value of the integer specified in such an I-field be j. Then the j'th molecular orbital in the reordered list will be of symmetry specified by the first data field of the current line. The following:
SIGMAG 1 2 3 4 5 6 7 8
when read as a 'molecular orbital ordering line' causes the first eight molecular orbitals of the reordered list to correspond to the first eight molecular orbitals of SIGMAG symmetry found in the input list of molecular orbitals. The 'molecular orbital ordering list' may be curtailed if sequences of consecutive integers appear, by means of the character string TO. Thus the above line is equivalent to:
SIGMAG 1 TO 8
Example :
Consider the hydrogen molecule, the basis set, as defined at integral evaluation time, being in the order:
H1(1s),H2(1s),H1(1s'),H2(1s')
Using the CTRANS directive, the user can symmetry adapt the basis functions:
CTRANS
1 1 2
1 2
1 3 2
1 4
1 1 2
-1 2
1 3 2
-1 4
END
To study the configuration (1-sigma-g)2(1-sigma-u)2 state of the H2 molecule, with a charge of -2, using the basis set defined above with the trial molecular orbitals generated through the START directive. The required configuration, as stated above, may be obtained using the data of the form:
ORBSYM
SIGMAG 1 2
SIGMAU 3 4
END
SIGMAG 1
SIGMAU 2
END
If the order of the trial molecular orbitals was 1-sigma-g, 2-sigma-g, 1-sigma-u, 2-sigma-u. After obeying the ORBSYM directive, the order would be 1-sigma-g, 1-sigma-u, 2-sigma-g, 2-sigma-u. Thus the first molecular orbitals found of sigma-g and sigma-u symmetry are placed in positions 1 and 2 in the re-ordered list. The remaining two molecular orbitals, which are not explicitly mentioned in the ORBSYM directive, are appended to positions 3 and 4, in the same order as they are found in the original configuration.
The ALIGN directive has been coded to allow the user to specify the alignment of degenerate sets of orbitals (usually occurs in GRHF module). The user supplies symmetry designation tags for the LCBF (see CTRANS directive). The program analyses each trial molecular orbital by means of a population analysis, and assigns symmetry designation tags to the molecular orbitals.
Suppose the user has labelled LCBF with the tags SIGMA, PIX and PIY. A given molecular orbital is found on Mulliken analysis to consist of 1% SIGMA, 20% PIX and 79% PIY, it will be labelled PIY. The program then deletes all 'contaminating' components from the molecular orbitals. Thus, in the given example, the coefficients of PIX and SIGMA orbitals will be set to zero. The GRHF module will retain the 'purified' state of the orbitals after the an ALIGN directive has been processed, by preventing the mixing of molecular orbitals of different symmetry. All other modules may allow such mixing however, and cause the iterated molecular orbitals to become 'inpure'.
The ALIGN directive may not appear before the LOADQ directive, and LCBF given different symmetry tags must have zero overlap integrals.
The first line of the ALIGN directive, the directive initiator, consists of the character string ALIGN, in the first data field. Following the directive initiator should be the basis function symmetry declaration lines, at least one line per irreducible representation being required. The first data field of such a line is read in A-format, and should contain a user specified symmetry label. Subsequent data fields are read in I-format; let the value of the integer specified in such an I-format field be j. Then the j'th basis function will be assigned the symmetry tag found in the first data field of that line. For example :
SIGMA 1 2 3 4 5 8 9 10 11
comprises a valid basis function symmetry declaration line, and assigns the tag SIGMA to LCBF 1 to 5 and 8 to 11. If sequences of consecutive integers arise, the line may be abbreviated using the character string TO. Thus :
SIGMA 1 TO 5 8 TO 11
is equivalent to the previous example. Any basis function not given a symmetry tag by the user will be assigned the tag '*' by the system. When all the basis function symmetry declaration lines have been read in, the directive is terminated by a data line containing the character string END. It is not possible to use the symmetry tag END. The ALIGN directive issues an ORTHOG directive in SCHMIDT mode.
Example :
A diatomic molecule has been oriented down the z-axis. The basis set, in the order, is, 1sa, 2sa, 3sa, 1sb, 2sb, 3sb, 2pxa, 2pya, 2pza, 2pxb, 2pyb, 2pzb. It is desired to align a trial set of molecular orbitals such that they are composed of only 'sigma', 'pi-x' or 'pi-y' components. The appropriate ALIGN directive would look like:
ALIGN
SIGMA 1 TO 5 8 11 14
PIX 6 9 12
PIY 7 10 13
END
an alternative form (allowing 'sigma' orbitals to be tagged '*') is:
ALIGN
PIX 6 9 12
PIY 7 10 13
END
This directive is used primarily to permit the SCF program to be capable of reconstruction of the eigen vectors as stored on the DUMP FILE from a file (card format). The MOPOPS directive will normally be followed by the STOP directive, when the latter will cause the eigen vectors restored by means of a LOADQ directive plus the molecular orbital populations defined by the MOPOPS directive to be written to a specified section of the DUMP FILE. The first line of the MOPOPS directive is read to variables TEXT,FORM using format (2A).
TEXT should be set to the character string MOPOPS.
FORM should be set to a valid FORTRAN format to be used for reading in the molecular orbitals occupation numbers. This parameter may be omitted, when the default 5F.16.9 will be chosen.
Subsequent data lines are read using the code:
READ(5,FORM)(POP(I),I=1,NEWBAS)
where FORM denotes the format specified by the FORM parameter (or 5F16.9 if FORM is omitted).
This directive is used primarily to permit the SCF program to be capable to reconstruct the eigen vectors as stored on the DUMP FILE, this directive is similar to the MOPOPS directive. The EVALUES directive allows eigen values (orbital energies) to be input from a file (card image) and subsequently to be written to the DUMP FILE. The specified section being designated by the STOP directive, which normally follows the EVALUES directive. The first line of the EVALUES directive is read to variables TEXT,FORM using format (2A).
TEXT should be set to the character string EVALUES.
FORM should be set to a valid FORTRAN format to be used for reading the eigen values. This parameter may be omitted, when the default 5F16.9 will be chosen.
Subsequent data lines are read using the code:
READ(5,FORM)(EIG(I),I=1,NEWBAS)
Note that the parameters defined by the MOPOPS or EVALUES directives take no active part in the SCF iteration process (indeed they will be overwritten during the iteration phase) such that these directives are solely for the purpose of reconstruction of the DUMP FILE.
This directive consists of a single data line read to variables TEXT, ISECT using format (A,I).
TEXT should be set to the character string ENTER.
ISECT is an integer specifying the section number of the DUMP FILE where the iterated eigen vectors are to be placed.
The actions of the ENTER directive are:
The ENTER directive must be presented last, any directives presented after the ENTER directive will be ignored. The ENTER directive is illegal if a LOADQ directive has not been found in the preceding data.
Example :
ENTER 10
This directive may be used as an alternative to the ENTER directive, and consists of one data line read to variables TEXT,ISECT using format (A,I).
TEXT should be set to the character string STOP.
ISECT an integer specifying the section number of the DUMP FILE where the trial vectors are to be written.
The effect of the STOP directive is as for the ENTER directive described previously, except that the SCF iterations are by-passed, execution being terminated as soon as the trial vectors have been written to the appropriate section of the DUMP FILE. The major use of the STOP directive is where the user wishes either:
The following notes are based on a paper by V. R. Saunders and I. H. Hillier [5]. The SCF calculation is performed in a basis of NEWBAS linearly independant LCBF. For a 2*NCORE electron system, one constructs a set of NCORE doubly occupied molecular orbitals (DOMOS) and NVIRT = NEWBAS - NCORE virtual molecular orbitals (VMOS), such that the complete set of molecular orbitals is orthonormal. Let 0, X1 and X2 denote row vectors of the LCBF, DOMOS and VMOS respectively.
(X1:X2) = 0(Q1:Q2) = 0Q 4.1
where a column of Q contains the atomic orbital coefficients of a given molecular orbital. Previously referred to Q as the 'eigen vector array'. The total wavefunction is invariant to mixing between DOMOS, so that mixing between DOMOS and VMOS, need only be of concern. Small variations may be written:
X1' = (X1:X2) !I1! 4.2
!..!
!D !
where I1 denotes an identity matrix of order NCORE, and D is the NVIRT*NCORE mixing matrix, whose elements are arbitrary but presumed small. The perturbed DOMOS are not orthonormal in second or higher order. Fortunately, the expansion of the energy is correct to first order only, which is given by:
virt occ
E = E0 + 4 [ [ Dki Hki + higher terms 4.3
k i
where E0 is the energy of the unperturbed wavefunction, and Hki denotes the matrix element connecting the k'th VMO with the i'th DOMO over the Roothaan-Hartree- Fock operator. From the form of equation 4.1 that the conditions of self consistency are satisfied when all Hki are zero.
! dE ! ------ D = 0 = 4Hki 4.4 !dDki!
It can be seen from 4.3 that any method capable of generating values for the elements of D such that
From the trial set of DOMOS, the Roothaan-Hartree-Fock operator is constructed in the basis set matrix representation and diagonalised. The resulting eigen vectors, after ordering, usually based on the aufbau principle, define improved molecular orbitals as linear combinations of the basis functions. A minor variant of the above, yielding equivalent results, is to diagonalise the Roothaan-Hartree-Fock operator constructed in the basis of the trial molecular orbitals. The resulting eigen vectors will define improved molecular orbitals as linear combinations of the trial molecular orbitals. The Roothaan- Hartree- Fock matrix in the trial molecular orbital basis can be partition as follows:
!Hdd:Hdv!
H = !-------! 4.5
!Hvd:Hvv!
where the partitions are defined as follows:
Hdd is the Fock matrix in the basis of the DOMOS
Hvv is the Fock matrix in the basis of the VMOS
Hdv (and Hvd) represents the Fock matrix connecting DOMOS with VMOS
The closed shell SCF program constructs a modified Fock matrix, Hmod, in the trial molecular orbital basis,
!Hdd : lHdv !
Hmod = !----:-------! 4.6
!lHvd:Hvv+aI2!
Hmod differs from H (4.5) because elements of the off-diagonal DOMO-VMO block have been multiplied by the damp factor (l), whilst the level shifter (a) has been added to each diagonal element of Hvv. I2 denotes the identity matrix of order NVIRT. The procedure adopted is to diagonalise the modified Fock matrix, and order the eigen vectors according to the aufbau principle based on the resulting eigen values. The first NCORE vectors define iterated DOMOS as linear combinations of the trial molecular orbitals.
In the following it is assumed, without loss of generality, that the canonicalisation of the trial DOMOS and VMOS is such that Hdd and Hvv are rendered diagonal. An analysis of the diagonalisation of Hmod by first order Rayleigh-Schrodinger theory yields the result:
Dki = lHki / (Hii - Hkk - a) 4.7
where the matrix D was defined by 4.2. The damp factor (l) is normally set to the value unity, and therefore concentrate on the effect of the level shifter (a) on convergence, assuming a unit damp factor. If a is chosen positive and sufficiently large, then:
In summary, a sufficiently level shifted procedure will give convergence to a stationary point on the energy surface, given any trial set of molecular orbitals.
(k+1) (k) (k) (1) (1) (2) (3) (k)
Q = Q T = Q T T T ........T 4.8
4.8 is not numerically satisfactory, since the build up in round-off error
is accumlative over k cycles. Round-off error can occur both in the
indicated matrix multiplications, and in the diagonalisation used to
generate the T matrices. To stabilise this method, the Schmit
orthonormalisation procedure to refresh the orthonormality of the columns of
Q (over the matrix of the overlap integrals of the LCBF) is used. This
orthonormalisation phase has the side benefit of removing errors due to the
user suppling a non-orthonormal set of initial vectors.The directives for the 'half closed' shell restricted Hartree-Fock module (RHF) it is suggested that they are presented in the following order:
See the TITLE directive of closed shell program.
See the MAXCYC directive of closed shell program.
See the MFILE directive of closed shell program.
See the SIZE directive of closed shell program.
See the CTRANS directive of closed shell program.
See the FPRINT directive of closed shell program.
See the PUNCH directive of closed shell program.
See the FORMAT directive of closed shell program.
See the OUTPUT directive of closed shell program.
See the IPRINT directive of closed shell program.
See the ACCURACY directive of closed shell program.
See the TIME directive of closed shell program.
This directive consists of one data line read to variables TEXT,DS1, SV1,IBRK,DS2,SV2 using format (A,2F,I,2F).
TEXT should be set to the character string LEVEL.
DS1,SV1 are the Doubly Occupied Molecular Orbital (DOMO) - Singly Occupied Molecular Orbital (SOMO) and SOMO - Virtual Molecular Orbital (VMO) level shifters, will be set to the values DS1 and SV1 respectively.
IBRK is the iterative cycle number to which the DS1 and SV1 level shifters hold to.
DS2,SV2 is the DOMO-SOMO and SOMO-VMO level shifters after cycle IBRK.
An alternative form of the LEVEL directive is allowed, where the variables are read to TEXT,DS1,SV1 using format (A,2F). In this form, TEXT,DS1 and SV1 have their usual meanings, whilst the program will set IBRK=999 and DS2,SV2 to zero. The LEVEL directive may be omitted, when the following default settings will apply:
DS1,SV1,DS2,SV2 = 0.0 IBRK = 999
The following notes may prove helpful to the user if level shifters are used to aid convergence.
L(k) = E(k) for DOMOS
L(k) = E(k) +a for SOMOS
L(k) = E(k) +a+b for VMOS
where 'a' and 'b' denote the DOMO-SOMO and SOMO-VMO level shifters
respectively. It is often found that a reasonable rate of convergence is
obtained when the level shifters are chosen such that:
L(lowest energy VMO)-L(highest energy SOMO) = 0.5 Hartree
L(lowest energy SOMO)-L(highest energy DOMO) = 0.4 Hartree
The above equations may be re-written:
E(lowest energy VMO)-E(highest energy SOMO)+b = 0.5 Hartree
E(lowest energy SOMO)-E(highest energy DOMO)+a = 0.4 Hartree
If a calculation is found to diverge and that the IPRINT directive was used
in the divergent run. Then the non-level shifted orbital energies, E(k),
will be available, and may be used in conjunction with the above formulae
to determine reasonable values for the level shifters.
LEVEL 0 0.3
Example 1:
Consider the following 'non-level shifted' orbital energy spectrum for a triplet state:
... -0.7 -0.5 -0.3 -0.1 0.1 0.2 ...
DOMOS SOMOS VMOS
To increase the SOMO-DOMO energy gap to 0.4 Hartree it is required that the first level shifter to be 0.2 and to increase the VMO-SOMO energy gap to 0.5 Hartree, it would require a value for the second level shifter of 0.3 Hartree. Thus the LEVEL directive would be:
LEVEL 0.2 0.3
This directive may be used to set the DOMO-SOMO 'damp factor' (referred as lds), and consists of one data line, read to variables TEXT,D12 using format (A,F).
TEXT should be set to the character string D12.
D12 should be set to the required value of the lds parameter.
The D12 directive is normally omitted, when the value of unity will be assigned to this parameter. Users are directed to reference [7] to where the D12 directives proves useful.
This directive may be used to set the DOMO-VMO 'damp factor' (referred as ldv), and consists of one data line, read to variables TEXT,D13 using format (A,F).
TEXT should be set to the character string D13.
D13 should be set to the required value of the ldv parameter.
The D13 directive is normally omitted, when the damp factor is assigned the default value of unity.
This directive may be used to set the SOMO-VMO 'damp factor' (referred as lsv), and consists of one data line, read to variables TEXT,D23 using format (A,F).
TEXT should be set to the character string D23.
D23 should be set to the required value of the lsv parameter.
The D23 directive is normally omitted, then the damp factor is assigned the value of unity.
This directive consists of one data line, read to variables TEXT,HD using format (2A).
TEXT should be set to the character string DOMOS.
HD should set to one of the character strings H1,H2,H3,H4 or H5.
The Fock operator used for the canonicalisation of the DOMOS, will be chosen according to the specification of the HD parameter from H1 to H5 as detailed in table 5.1. The DOMOS directive may be omitted, when H1 will be used to canonicalise the DOMOS.
Example :
DOMOS H2
This directive consists of a single data line, read to variables TEXT, HS using format (2A).
TEXT should be set to the character string SOMOS.
HS should be set to one of the character strings H1,H2,H3,H4 or H5.
The Fock operator used for the canonicalisations of the SOMOS, is chosen according to the specification of the HS parameter, from H1 to H5, as outlined in table 5.1. The SOMOS directive may be omitted, when H1 will be used to canonicalise the SOMOS.
Example :
SOMOS H2
This directive consists of a single data line, read to variables TEXT, HV using format (2A).
TEXT should be set to the character string VMOS.
HV should be set to one of the character strings H1,H2,H3,H4 or H5.
The Fock operator used for the canonicalisation of the VMOS, is chosen to the specification of the HV parameter, from H1 to H5, as outlined in table 5.1.
Example :
VMOS H2
Table 5.1
_________
Definition of Fock Operators for the RHF Module
_______________________________________________
Fock Operator Cn
_____________ __
H1 0
H2 -1
H3 +1
H4 -2
H5 +2
This directive consists of a single data line with the character string LOCK in the first data field. If used, the LOCK directive causes iterated molecular orbitals to be selected on the principle of maximum overlap with trial molecular orbitals.
The LOCK directive may be omitted, when ordering is based on the 'aufbau' principle, so that the molecular orbitals are ordered according to increasing eigen value. The directive will be found useful for the study of excited states. The user should note that in open shell calculations, there are often a number of states close to the ground state energy. The user should not assume that because the self consistent solution obeys the 'aufbau' principle, it corresponds to the ground state. Conversely, a number of cases are known where the ground state does not obey the 'aufbau' principle, particularly where other than default canonicalisations have been specified (by means of the DOMOS,SOMOS or VMOS directives). The user should obtain convergence to one stationary point on the energy surface. If the resulting configuration is not the one required, convergence to another configuration may be achieved by restart runs employing SWAP and LOCK directives.
It is possible to run a closed shell calculation through the RHF module (by specifying NOPEN as zero in the first data line). The only reason for doing so is to make use of the LOCK option, the latter being not available with the SCF module.
One of the LOADQ directives VECTORS,START,ALPHAS,RESTORE,GETQ,EXTRA or COMBINE must be presented so that trial eigen vectors may be obtained. The first NCORE orbitals will be doubly occupied, the next NOPEN will be singly occupied, the remainder will be virtual.
See the SDIAG directive of closed shell program.
See the ORTHOG directive of closed shell program.
See the NOORTH directive of closed shell program.
See the SDWAP directive of closed shell program.
See the ORBSYM directive of closed shell program.
See the ALIGN directive of closed shell program.
See the MOPOPS directive of closed shell program.
See the EVALUES directive of closed shell program.
See the SENTER directive of closed shell program.
See the STOP directive of closed shell program.
The following notes are based on the paper by M. F. Guest and V. R. Saunders [8]. The purpose of the RHF program is to make stationary (usually to minimise) the energy of a single determinatal wavefunction contructed from NCORE orthonormal doubly occupied molecular obitals (DOMOS) and NOPEN orthonormal singly occupied molecular orbitals (SOMOS), latter being of common spin factor. From NEWBAS linearly independent basis functions one can construct NVIRT = NEWBAS-NCORE-NOPEN virtual molecular orbitals (VMOS) such that all molecular orbitals comprise an orthonormal set. Let 0, X1, X2 and X3 denote row vectors of the basis functions DOMOS, SOMOS and VMOS respectively. Then:
(X1:X2:X3) = 0(Q1:Q2:Q3) = 0Q 5.1
where each column of Q contains the basis function coefficients of a given molecular orbital. The density matrices R1 and R2 are defined as:
R1 = Q1Q1+ 5.2 R2 = Q2Q2+ 5.3
The electronic energy expression for the half-closed shell case is given by:
E = 2trR1F + trR2F + 2trR1J[R1] - trR1K[R1] + 1/2trR2J[R2]
- 1/2trK2[R2] + 2trR2J[R1] - trR2[R1] 5.4
where F is the matrix representative of the one-electron Hamiltonian operator; J[R] and K[R] denote Coulomb and exchange matrices respectively for a given density matrix, R. The functional notation is to be taken as meaning that the matrix elements of the J and K are functions of the matrix elements of R, as follows:
Jpg = [ Rrs(pq!rs) 5.5
rs
Kpq = [ Rrs(pr!qs) 5.6
rs
where
(pq!rs) = // 0p(1)0q(1) 1/r12 0r(2)0s(2) dt1 dt2 5.7
Now define:
R = 2R1 + R2 5.8 P = 1/2R2 5.9
and
G = J[R] - 1/2K[R] 5.10 Z = K[P] 5.11
It is convenient to make the following definition of a Fock operator:
Hn = F + G + cnZ 5.12
where the values for cn of interest are displayed in table 5.1. The energy expression (5.4) may be written in the computationally convenient form:
E = 1/2trR(F+H1) - trPZ 5.13
Physically significant variations in the molecular orbitals may be classified as follows:
X1' = (X1:X3) !I1! 5.14
!--!
!A !
where I1 is an identity matrix of order NCORE, and A is a NVIRT*NCORE
matrix of mixing coefficients.
X2' = (X2:X3) !I2! 5.15
!--!
!B !
where I2 is an identity matrix of order NOPEN, and B is a NVIRT*NOPEN
matrix of mixing coefficients.
(X1:X2)' = (X1:X2) !I1:-C+! 5.16
!--:---!
!C : I2!
where C is a NOPEN*NCORE matrix of mixing coefficients.Combining 5.14, 5.15 and 5.16,
(X1:X2)' = (X1:X2:X3) !I1:-C+! 5.17
!--:---!
!C : I2!
!--:---!
!A : B !
Notice that self mixing between DOMOS or SOMOS has no effect on the wavefunction of the energy. The perturbed molecular orbitals are not orthonormal in second or higher order. Fortunately, the energy expression is only required to the first order in the elements of the matrices A, B and C, given by:
NVIRT NCORE NVIRT NOPEN
E = E0 + 4 [ [ Aki H1ki + 2 [ [ Bkj H2kj
k i k j
NCORE NOPEN
= 2 [ [ Cji H3ji + higher terms 5.18
i j
H1ki denotes the matrix elements connecting the k'th VMO with the
i'th DOMO over the Fock operator H1 (see table 5.1).
H2kj denotes the matrix elements connecting the k'th VMO with the
j'th SOMO over the Fock operator H2 (see table 5.1).
H3ji denotes the matrix elements connecting the j'th SOMO with the
i'th DOMO over the Fock operator H3 (see table 5.1).
It is clear that the conditions for a stationary energy are satisfied when all matrix elements (H1ki, H2kj and H3ji) appearing in 5.18 are zero. Thus:
(dE ) -----Aki = 0 = 4H1ki 5.19 (dAki) (dE ) ------ Bkj = 0 = 2H2kj 5.20 (dBkj) (dE ) ------ Cji = 0 = 2H3ji 5.21 (dCji)
A suitable Fock operator id defined in the basis of the trial molecular orbitals (normally the Fock operator is defined in the basis function representation), and this Fock operator is diagonalised. The resulting eigen vectors, after suitable ordering, define iterated molecular orbitals as linear combinations of the trial molecular orbitals.
Let Q(k) denote the eigen vector array defining trial molecular orbitals as linear combinations of the LCBF at the k'th iterative cycle, and let T(k) denote eigen vectors resulting from the diagonalisation of the Fock operator in the iterated molecular orbital basis, hence defining iterated molecular orbitals as linear combinations of the trial molecular orbitals.
The general form of the Fock operator is given by:
! (HD)dd:lds(H3)ds :ldv(H1)dv !
----------------------------------------!
!lds(H3)sd:(HS)ss + aI2:lsv(H2)sv !
----------------------------------------!
!ldv(H1)vd:lsv(H2)vs :(HV)vv + (a+b)I3!
(HD)dd, (HS)ss and (HV)vv : each denote a block of one of the operators H1 to H5 in the basis of the DOMOS, SOMOS and VMOS respectively.
lds, ldv and lsv : denote DOMO-SOMO, DOMO-VMO and SOMO-VMO damp factors respectively.
a and b : denote the DOMO-SOMO and SOMO-VMO level shifters respectively.
(H3)ds : denotes a block of H3 connecting all DOMOS with all SOMOS. (H3)sd is the transpose.
(H1)dv : denotes a block of H1 connecting all DOMOS with all VMOS. (H1)vd is the transpose.
(H2)sv : denotes a block of H2 connecting all SMOS with all VMOS. (H2)vs is the transpose.
It will be readily seen that if the above Fock operator is diagonal, the conditions for a stationary energy will be satisfied, so that the total wavefunction will be self-consistent. The procedure adopted is to diagonalise the Fock operator, and unless the LOCK directive is presented, to order the resulting eigen vector columns according to the level shifted eigen values. If the LOCK directive was used, then ordering is based on the principle of maximum overlap of an iterated eigen vector with the trial eigen vector. Whichever ordering scheme is adopted, the first NCORE and the next NOPEN columns define iterated DOMOS and SOMOS as linear combinations of the trial molecular orbitals.
Analyse of the diagonalisation of the Fock operator using first order perturbation theory, yields the following results:
Aki = ldv H1ki / (HDii - HVkk - a - b) 5.23 Bkj = lsv H2kj / (HSjj - HVkk - b) 5.24 Cji = lds H3ji / (HDii - HSjj - a) 5.25
The following discusion is based on the assumption that the usual value of unity has been chosen for the l-parameters (ldv, lds and lsv). If a and b are chosen positive and with a sufficiently large value, then:
A sufficiently level shifted procedure will give convergence to a stationary point on the energy surface given any trial set of molecular orbitals, and will always procede down the energy surface.
The level shifted eigen values of the Fock operator are reduced by the constants (a+b) and a for the VMOS and SOMOS respectively, so as to produce non-level shifted eiegn values on the printed output.
At convergence, the canonicalisation of the molecular orbitals is such that:
Q1+(HD)Q1 = Ad 5.26 Q2+(HS)Q2 = As 5.27 Q3+(HV)Q3 = Av 5.28
where Ad, As and Av are diagonal matrices defining the non-level shifted orbital energies of the DOMOS, SOMOS and VMOS respectively. The purpose of the following is to give some idea as to why other than the default canonicalisation scheme (which normally gives a good convergence rate) might be used.
The Roothann Fock operator which appears in Roothaan's one-electron Hamiltonian formulation [4], may be written in the basis of the trial molecular orbitals, and takes the form:
!(H4)dd:(H3)ds:(H1)dv!
!--------------------!
!(H3)sd:(H1)ss:(H2)sv!
!--------------------!
!(H1)vd:(H2)vs:(H5)vv!
it will be readily seen that the Roothaan single Fock operator is a special case of the ATMOL Fock operator, and the program may be forced to perform iterations in precisely the manner prescribed by Roothaan's single Fock operator method if the directives:
DOMOS H4
SOMOS H1
VMOS H5
are presented in the user's input data.
The Roothaan double Fock operator method cannot be shown to be a special case of the ATMOL procedure. However identical orbitals and orbital energies result from the present program if canonicalisations of the DOMOS and SOMOS are specified as H1 and H2 respectively.
The canonicalisations may be specified so that the orbital energies of the DOMOS and SOMOS are ionisation potential in the Koopmans' theorem sense. Assume that the common spin factor of the SOMOS is a, and consider the removal of an electron from the b-spin DOMO, giving rise to an ionic configuration of spin degeneracy one greater than the parent molecular configuration. It is possible to form NCORE such configurations, and would expect in general, the description of the ionic states would be improved by allowing configuartion interaction (CI) between these configurations. The ion Hamiltonian matrix element connecting two ionic configurations, Xp and Xq, may be expressed as:
Hpq = = - 5.29
where Wp and Wq denote the molecular orbitals which have under-gone electron removal to form the configurations, Xp and Xq. Note that if H3 had been chosen to canonicalise the DOMOS, then Hpq would be equal to zero, and no interaction between the above ionic configurations would occur. Moreover it can be easily shown that the energy required for the removal of an electron from b-spin DOMO is given by the negative of the orbital energy of the DOMO, calculated over H3.
Now consider the case of removal of an electron from an a-spin SOMO, giving rise to an ionic configuration of spin degeneracy one less than the parent molecular configuration. It is possible to form NOPEN such configurations, and again would expect that the description of the ionic states would be improved by allowing CI between the configurations. The Hamiltonian matrix connecting two configuartions, Xq and Xv is given by:
Hqv = = - 5.30
where Wq and Wv denote the molecular orbitals which have under-gone electron removal to form the configurations Xq and Xv. Thus, if H2 had been chosen to canonicalise the DOMOS, then Hqv would equal zero, and no interaction between the above configurations would occur. Moreover it can be shown that the energy required for the removal of an electron from a SOMO is given by the negative of the orbital energy of the SOMO, calculated over H2.
The canonicalisation of the SOMOS and VMOS may be chosen so as to produce electron affinity orbital energies. Consider the addition of an electron to a SOMO of b-spin, to produce an ionic configuration whose spin degeneracy is one less than that of the parent molecule. NOPEN such configurations may be produced, and in general would expect the description of the ion states to be improved by CI. The matrix element connecting two ionic configurations Xp and Xq is given by:
= 5.31
where Wp and Wq denote the SOMOS involved as electron acceptors in the configurations Xp and Xq respectively. Thus, if the SOMOS are canonicalised over H3, the ionic configurations formed by addition of an electron to a SOMO will be non-interacting. Furthermore, the energy released when an electron is added to the p'th SOMO (electron affinity) is given by the negative of the corresponding orbital energy computed over H3, ep. Thus:
Eion = Emol + ep 5.32
Consider the addition of an electron to an a-spin VMO, to produce an ionic configuration whose spin degeneracy is one greater than that of the parent molecule. NVIRT such configurations can be generated, and Hamiltonian matrix elements are given by:
= 5.33
where Wp and Wq denote VMOS involved as electron acceptors in the configurations Xp and Xq respectively. Thus, if the VMOS are canonical over H2, the above ionic configurations will be non-interacting. Further, the energy released when an electron is added to the p'th a-spin VMO is given by the negative of the corresponding orbital energy, computed over H2, ep. Thus:
Eion = Emol + ep 5.34
Table 5.2
_________
Useful Canonicalisations
Optimum Ionisation Electron To simulate To reproduce
Convergence Potentials Affinities Roothaan's results from
Single Fock Roothaan's Double
Matrix Method Fock Matrix Method
------------------------------------------------------------------------------
DOMOS H1 H3 Not Valid H4 H1
SOMOS H1 H2 H3 H1 H2
VMOS H1 Not Valid H2 H5 Not Valid
The aim of the symmetry equivalenced RHF module (codenamed SERHF) is to provide facilities for the minimisation of energy expressions which are of more general form than can be handled by the RHF module. This module can be used to symmetry equivalence orbitals, such that an average energy expression may be obtained. The directives of the SERHF module are now outlined, generally they may be presented in any order if not previously stated.
See the TITLE directive of closed shell program.
See the MAXCYC directive of closed shell program.
See the MFILE directive of closed shell program.
See the SIZE directive of closed shell program.
See the CTRANS directive of closed shell program.
See the FPRINT directive of closed shell program.
See the PUNCH directive of closed shell program.
See the FORMAT directive of closed shell program.
See the OUTPUT directive of closed shell program.
See the IPRINT directive of closed shell program.
See the ACCURACY directive of closed shell program.
See the TIME directive of closed shell program.
This directive consists of one data line read to variables TEXT,DP1,PV1, IBRK,DP2,PV2 using format (A,2F,I,2F).
TEXT should be set to the character string LEVEL.
DP1,PV1 are the Doubly Occupied Molecular Orbital (DOMO)-Partially Occupied Molecular Orbital (POMO) and POMO - Virtual Molecular Orbital (VMO) level shifters, will be set to the values DS1 and SV1 respectively.
IBRK is the iterative cycle number to which the DP1 and PV1 level shifters hold to.
DP2,PV2 is the DOMO-POMO and POMO-VMO level shifters after cycle IBRK.
In general terms, the partially occupied molecular orbitals (POMO) have a similar role within the SERHF module as the SOMOS have for the RHF module. Users may find it useful to read the description of the LEVEL directive for the RHF program to gain insight into the corresponding use within the SERHF program.
These directives are used to set the parameters ldp, ldv and lpv respectively referred to in section 6.xx. The general syntax may be deduced from the descriptions given for the corresponding directives of the RHF program (paragraphs 5.xx t0 5.yy)
This directive is used to define the energy expression to be minimised, and consists of one data line, read to variables TEXT,NELEC,X,Y using format (A,I,2F).
TEXT should be set to the character string ENERGY.
NELEC should be set to the total number of electrons in the open shell orbitals. Defined in a subsequent table is the parameter 'w', such that w = NELEC / NOPEN.
This directive consists of a single data line, read to variables TEXT, CD,DD using format (A,2F).
TEXT should be set to the character string DOMOS.
CD,DD are used to specify the values of the parameter 'c' and 'd' respectively (see 6.5), which enters into the canonicalising Fock matrix for the DOMOS, referred as HD in paragraph 6.xx.
The DOMOS directive may be omitted, when H1 (as defined in table 6.3) will be used to canonicalise the DOMOS.
This directive consists of one data line read to variables TEXT,CP,DP using format (A,2F).
TEXT should be set to the character string POMOS.
CP,CD are used to specify the values 'c' and 'd' respectively (see 6.5) which enter into the canonicalising Fock matrix for the POMOS, referred to as HP in paragraph 6.xx.
The POMOS directive may be omitted, when the POMOS will be canonicalised over H1, as defined in table 6.3.
This directive consists of a single data line, read to variables TEXT, CV,DV using format (A,2F).
TEXT should be set to the character string VMOS.
CV,CD are used to specify the values 'c' and 'd' respectively, (see 6.5) which enter into the canonicalised Fock matrix for the VMOS, referred to as HV in paragraph 6.xx.
The VMOS directive may be omitted, when the VMOS will be canonicalised over H1, as defined in table 6.3.
See the LOCK directive of open shell program.
One of the LOADQ directives (VECTORS,START,ALPHAS,RESTORE,GETQ,EXTRA or COMBINE) must be presented, to define trial eigen vectors. The first NCORE orbitals will be doubly occupied, the next NOPEN will be partially occupied, the remainder virtual.
See the SDIAG directive of closed shell program.
See the ORTHOG directive of closed shell program.
See the NOORTH directive of closed shell program.
See the SDWAP directive of closed shell program.
See the ORBSYM directive of closed shell program.
See the ALIGN directive of closed shell program.
See the MOPOPS directive of closed shell program.
See the EVALUES directive of closed shell program.
See the SENTER directive of closed shell program.
See the STOP directive of closed shell program.
The following notes are based on a paper by M.F Guest and V.R Saunders [8]. The purpose of the SERHF program is to make stationary (often to minimise) the energy of a wavefunction constructed from NCORE doubly occupied molecular orbitals (DOMOS) and NOPEN partially occupied molecular orbitals (POMOS). From NEWBAS linearly independent basis functions one can construct NVIRT = NEWBAS - NCORE - NOPEN virtual molecular orbitals (VMOS) such that all the molecular orbitals comprise an orthonormal set. Let 0, X1, X2, and X3 denote row vectors of the basis functions, DOMOS, POMOS and VMOS respectively. Then:
(X1:X2:X3) = 0(Q1:Q2:Q3) = 0Q 6.1
where each column of Q contains the basis function coefficients of a given molecular orbital. The density matrices are defined as:
R1 = Q1 Q1+ 6.2 R2 = Q2 Q2+ 6.3
The electronic energy expression to be minimised is assumed to be of the form:
E = 2trR1F + wtrR2F + 2trR1J[R1] -trR1K[R1]
+ 2wtrR2J[R1] - wtrR2K[R1] + xtrR2J[R2] - ytrR2K[r2] 6.4
where w, x and y are the parameters discussed in subparagraph 6.xx, F is the matrix representative of the one-electron Hamiltonian, whilst the construction of the Coulomb and exchange matrices, J[R] and K[R] respectively, have been outlined in subparagraph 5.33.1.
It is convenient to make the following definition of a class of Fock operators:
H = F + 2J1[R1] - K[R1] + cJ[R2] - d1/2K[R2] 6.5
In table 6.3 the c and d values are given (as functions of the energy expression parameters w, x and y) for those forms of H that the user may find useful.
Table 6.3
_________
Definition of Fock Operators for the SERHF Module
_________________________________________________
Fock Operator c d
_____________ _ _
H1 w w
H2 2x/w 4y/w
H3 2(w-x)/(2-w) 2(w-2y)/(2-w)
2 3 2 3
H4 (4w -w -4x)/w(2-w) (4w -w -8y)/w(2-w)
3 3
H5 (w -4xw+4x)/w(2-w) (w -8yw+8y)/w(2-w)
3 3
H6 (4x-w )/w(2-w) (8y-w )/w(2-w)
Physically significant variations in the molecular orbitals may be classified as follows:
Mixing of DOMOS with VMOS. Small such variations preserving the orthonormality of the DOMOS to first order, may be written:
W1' = (W1:W3) !I1!
!--! 6.6
!A !
where I1 is an identity matrix of order NCORE, and A is a NVIRT*NCORE matrix of mixing coefficients.
Mixing of POMOS with VMOS. Small such variations, preserving the orthonormality of the POMOS to first order, may be written:
W2' = (W2:W3) !I2!
!--! 6.7
!A !
where I2 is an identity matrix of order NOPEN, and B is a NVIRT*NOPEN matrix of mixing coefficients.
Mixing of DOMOS and POMOS. Small such variations, preserving the orthonormality of the DOMOS and POMOS to first order, may be written:
(W1:W2)' = (W1:W2) !I1:-C+!
!------! 6.8
!C : I2!
where C is a NOPEN*NCORE matrix of mixing coefficients.
Equations 6.6 to 6.8 may be combined to produce:
(W1:W2)' = (W1:W2:W3) !I1:-C+!
!------!
!C : I2! 6.9
!------!
!A : B !
Notice that self mixing between the DOMOS or POMOS has no effect on the wavefunction or energy. The perturbed molecular orbitals are not orthonormal beyond first order. The expansion of the energy expression is required to be correct only to the first order in the elements of the matrices A, B and C, given by:
NVIRT NCORE NVIRT NOPEN
E = E0 + 4 [ [ AkiH1ki + 2w [ [ BkjH2ki
k i k j
NCORE NOPEN
+ 2(2-w) [ [ CjiH3ji + higher terms 6.10
i j
H1ki denotes the matrix element connecting the k'th VMO with the i'th DOMO over the Fock operator H1 (see table 6.3).
H2kj denotes the matrix element connecting the k'th VMO with the j'th POMO over the Fock operator H2 (see table 6.3).
H3ji denotes the matrix element connecting the j'th POMO with the i'th DOMO over the Fock operator H3 (see table 6.3).
It is clear that the conditions for a stationary energy are satisfied when all the matrix elements (H1ki, H2kj and H3ji) appearing in 6.10 are zero.
The general strategy of the SERHF energy is as for the RHF program. The ATMOL SERHF Fock operator in the basis of the trial molecular orbital is given by:
! (HD)dd:lds(H3)dp : ldv(H1)dv !
!-----------------------------------!
!ldp(H3)pd:(HP)pp+aI2: lpv(H2)pv !
!-----------------------------------!
!ldv(H1)vd:lpv(H2)vp :(HV)vv+(a+b)I3!
where:
(HD)dd, (HP)pp and (HV)vv: denote a block of a Fock operator of the general form given in equation 6.5, in the basis of DOMOS, POMOS and VMOS repectively. The user may specify the form of HD, HP and HV as described in subparagraphs 6.xx to 6.xx.
ldp, ldv and lpv: denote the DOMO-POMO, DOMO-VMO and POMO-VMO damp factors, referred in subparagraph 6.xx.
a and b: denote the DOMO-POMO and POMO-VMO level shifters respectively.
(H3)dp: denotes a block of H3 (see table 6.3) connecting all DOMOS with all POMOS. (H3)pd is the transpose.
(H1)dv: denotes a block of H1 (see table 6.3) connecting all DOMOS with all VMOS. (H3)pd is the transpose.
(H2)pv: denotes a block of H2 (see table 6.3) connecting all POMOS with all VMOS. (H2)vp is the transpose.
The Fock operator which appears in the Roothaan's one-electron Hamiltonian [4], may be written in the basis of the trial molecular orbitals, and takes the form:
!(H4)dd:(H3)dp:(H1)dv!
!--------------------!
!(H3)pd:(H5)pp:(H2)pv!
!--------------------!
!(H1)vd!(H2)vp:(H6)vv!
It will be readily seen that the Roothaan single Fock operator is a special case of the ATMOL SERHF Fock operator, and the SERHF program may be forced to perform the SCF iterations in precisely the manner prescribed by Roothaan's single Fock operator, given that appropriate DOMOS, POMOS and VMOS directives are used.
The Roothaan double-Hamiltonian method cannot be shown to be a special case of the SERHF method. However, it can be shown that the orbital energies and canonicalisations of the converged DOMOS and POMOS may be rendered identical in the two procedures with the choice HD = H1, HP = H2.
It may be useful to note the following connections between the Roothaan f, a and b energy expression parameters, and the SERHF w, x and y parameters:
w = 2f 6.11 x = 2af2 6.12 y = bf2 6.13
The aim of the 'general' RHF module (code named GRHF) is to provide facilities for the minimisation of energy expressions of a more general form than can be handled by the RHF or SERHF modules. The user should note that the NCORE and NOPEN parameters of the first data line have no significance for the GRHF module, although they are submitted to certain standard checking procedures, and the users are thus advised to set both of these parameters to unity. The directives of the GRHF module following the first data line may be presented in any order, unless explicitly specified.
See the TITLE directive of closed shell program.
See the MAXCYC directive of closed shell program.
See the MFILE directive of closed shell program.
See the SIZE directive of closed shell program.
See the CTRANS directive of closed shell program.
See the FPRINT directive of closed shell program.
See the PUNCH directive of closed shell program.
See the FORMAT directive of closed shell program.
See the OUTPUT directive of closed shell program.
See the IPRINT directive of closed shell program.
See the ACCURACY directive of closed shell program.
See the TIME directive of closed shell program.
This directive, which is comparitively complex in structure, defines data which plays an equivalent role in the GRHF module to that supplied by the LEVEL,D12,D13,D23,DOMOS,POMOS,VMOS and ENERGY directives of the SERHF module. The SHELLS directive has the additional responsibility of defining the molecular orbital 'shell' structure.
The first data line should contain the character string SHELLS in the first data field. Subsequent data lines define the molecular orbital shell structure, one data line per shell, each line being read in I-free format, the integers specifying the numbering of the molecular orbitals to appear in a given shell. For example, if the following line appeared:
1 3 5
as the first 'shell definition line'. Then molecular orbitals 1,3 and 5 would be deemed to constitute shell 1. Similarly, if the data line:
10 11
appeared as the second 'shell definition line', then molecular orbitals 10 and 11 would be deemed to constitute shell 2. The use of the character string TO may be used to simplify the preparation of 'shell definition lines'. For example, the following data lines are equivalent in the context:
1 2 3 4 5 7
or
1 TO 5 7
or
7 1 TO 5
When all the 'shell definition lines' have been presented, the sequence should be terminated by a line containing the character string END in the first data line. A maximum of ten shells can be handled by the GRHF module.
The data lines following the 'shell definition lines' are used to define the 1-electron energy expression parameters (W parameters, see equation 3). The data is read in free F-format to an array, (W(I),I=1,NSHELLS), where NSHELLS equals the number of occupied shells. As many data lines as required may be used, so that sequences:
2.0 1.0 0.333333
and
2.0
1.0
0.333333
are equivalent.
Immediately following the 1-electron energy expression definition lines there follows a sequence of lines read to variables TEXT,M,N,VAL using format (A,2I,F). The purpose of these lines is to define coulomb and exchange energy expression parameters, damp factors and level shifters and canonicalisation conditions.
The role of each data line is dependent upon the character string found in TEXT, as follows:
TEXT = JE : Such lines define coulombic energy expression parameters.
The program will arrange that x = x , so that there is
NM MN
no need for the user to specify 'x' directly.
If the coulomb energy expression parameter for a given pair of shells is not user defined by means of a JE line, it will be given the default value:
x = (W W ) / 2 = x
MN M N NM
Note that the maximum value for M and N is NSHELL.
TEXT = KE: Such lines define the exchange expression parameters. The
program will arrange that y = y .
NM MN
If the energy energy expression parameter is not user defined by means of a KE line, it will be given the default value:
y = -(W W ) / 4
MN M N
Note that the maximum value for M and N is NSHELL.
TEXT = DAMP: Such lines define tha damp factor for a given pair of shells. If the damp factor for a given pair of shells is not user defined, it will be assigned the default value unity by the program.
Because mixing between occupied and virtual shells is permitted in the minimisation process, the maximum value of M and N is given by NSHELL+1.
TEXT = SHIFT: Such lines define the level shifter for a given pair of shells. If the level shifter for a given pair of shells is not user defined, it will be assigned the default value zero by the program.
The maximum value of M and N is given by NSHELL+1.
TEXT = JC: Such lines define coulomb canonicalisation factors. If the user does not define 'c' (see equation 4) explicitly, by means of the JC line, it will be assigned the value W by the program.
The maximum value of M is given by NSHELL+1 (the program does canonicalise VMOS), the maximum value of N being NSHELL (because only the occupied shells have their coulomb operators calculated). Note that the program does not assume or require the matrix of 'c' parameters to be symmetric.
TEXT = KC: Such lines define exchange canonicalisations factors. If the user does not define 'd' (see equation 4) explicitly, by the means of a KC line, it will be assigned the value -W/2 by the program.
The maximum value of M is given by NSHELL+1, the maximum value of N is NSHELL. The program does not assume or require that the matrix of 'd' parameters is symmetric.
When all the JE,KE,DAMP,SHIFT,JC and KC lines have been presented (in any order), the SHELLS directive is terminated by a line containing the character string END in the first data field.
One of the LOADQ directives (VECTORS,START,ALPHAS,RESTORE,GETQ,EXTRA or COMBINE) must be presented, to define trial eigen vectors.
See the SDIAG directive of closed shell program.
See the ORTHOG directive of closed shell program.
See the NOORTH directive of closed shell program.
See the SDWAP directive of closed shell program.
See the ORBSYM directive of closed shell program.
See the ALIGN directive of closed shell program.
See the MOPOPS directive of closed shell program.
See the EVALUES directive of closed shell program.
See the SENTER directive of closed shell program.
See the STOP directive of closed shell program.
The GRHF program assumes that the total wavefunction is constructed from NACT occupied molecular orbitals which are orthonormal. It further assumes that the user has assigned the occupied molecular orbitals to 'shells', where a 'shell' consists of a set of one or more molecular orbitals. Note that a given molecular orbital may be assigned to only one shell. In the following notes the symbols M and N denote shells, and the symbols m and n denote molecular orbitals within shells M and N respectively. From NEWBAS linearly independent basis functions the program is able to construct NEWBAS - NACT VMOS. Given that NSHELLS shells have been defined as occupied, then the VMOS will comprise the NSHELL+1'th shell, with occupation numbers of zero. The VMOS will be so constructed that the total set of NEWBAS molecular orbitals are orthonormal. The following parameters are defined as:
Fmm The expectation value of the m'th molecular orbital
over the unscreened 1-particle Hamiltonian operator.
FM = [ Fmm The sum of the 1-particle expectation values for all
molecular orbitals in shell M.
Jmn = [mm!nn] The Columb integral between molecular orbitals m
and n.
JMN = [ [ Jmn
Kmn = [mn!mn] The exchange integral between molecular orbitals m
and n.
KMN = [ [ Kmn
Energy expressions which can be handled by the GRHF program are of the form:
NSHELL NSHELL NSHELL NSHELL NSHELL
E = [ WM FM + [ [ xMN JMN + [ [ yMN KMN 7.1
M M N M N
Where WM are the one electron energy expression parameters.
xMN are the Coulomb energy expression parameters (note that xMN must equal xNM).
yMN are the exchange energy expression parameters (note that yMN must equal yNM).
The SHELLS directive is used to define the above energy expression parameters. Note the form of equation 7.1, that the energy is invariant to mixing of the molecular orbitals within a shell.
Detailed below are simple examples of energy expressions which can be handled by the GRHF program. More detailed examples of particular energy expressions of selected states are given at the end of this section.
Example 1
Consider the case of a half-closed shell triplet, the wavefunction may be written as:
- -
!a a b b c d!
In this case the orbitals a and b comprise one shell (the closed shell) and orbitals c and d comprise another (the open shell). The energy expression for such a state (see equation 5.4) may be rewritten in the form:
E = 2F1 + F2 + 2J11 +J12 + J21 + 1/2J22 - K11 - 1/2K12 - 1/2K21 - 1/2K22
The SHELLS directive to describe such a case would look like:
SHELLS
1 2
3 4
END
2.0 1.0
JE 1 1 2.0
JE 1 2 1.0
JE 2 2 0.5
KE 1 1 -1.0
KE 1 2 -0.5
KE 2 2 -0.5
END
However, remember that the program sets Coulombic and exchange energy expressions to default values according to the rule:
xMN = (WMWN) / 2 7.2 yMN = -(WMWN) / 2 7.3
The only parameter which does not comply with the above default settings is that for exchange self energy of shell 2. Thus the above SHELLS directive may be simplified to the form:
SHELLS
1 2
3 4
END
2.0 1.0
KE 2 2 -0.5
END
Example 2:
Consider the following singlet wavefunction:
_ _ _ _ _ _
(!a a b b c d! + !a a b b d c!) / (2)**0.5
In this case, orbitals a and b comprise one shell (the closed shell), orbital c comprises the second shell and orbital d comprises the third shell. The energy expression may be written:
E = 2F1 + F2 + F3 + 2J11 + J12 + J21 + J13 + J31 + 1/2J22
+ 1/2J23 + 1/2J32 + 1/2J33 - K11 - 1/2K12 - 1/2K21 - 1/2K13
- 1/2K31 - 1/2K22 + 1/2K23 + 1/2K32 - 1/2K33
The appropriate SHELLS directive to define the above case would look like:
SHELLS
1 2
3
4
END
2.0 1.0 1.0
KE 2 2 -0.5
KE 2 3 0.5
KE 3 3 -0.5
END
The only two-electron energy expression parameters which differ from the defaults are those for the exchange self energies of shells 2 and 3, and also the exchange energy between shells 2 and 3.
Example 3:
Consider the two term configuration interation (CI) function:
_ _ _ _ _ _
p!a a b b c c! + q!a a b b d d!
where p2 + q2 = 1
Normally p and q will have been determined from a CI calculation. Let molecular orbitals a and b constitute shell 1, orbital c constitute shell 2 and orbital d constitute shell 3. The one-electron energy expression parameters are given as:
W1 = 2 : W2 = 2p2 : W3 = 2q2
The Coulomb energy expression parameter matrix written:
SHELLS 1 2 3 ------------------------------------------ ! 1 ! 2 ! 2p2 ! 2q2 ! ------------------------------------------ ! 2 ! ! 2p2 ! 0 ! ------------------------------------------ ! 3 ! ! ! 2q2 ! ------------------------------------------
where the lower triangle of the matrix has been omitted (the matix is symmetric).
The exchange energy expression parameter matrix may be written as:
SHELLS 1 2 3 ------------------------------------------ ! 1 ! -1 ! -p2 ! -q2 ! ------------------------------------------ ! 2 ! ! -p2 ! pq ! ------------------------------------------ ! 3 ! ! ! -q2 ! ------------------------------------------
again the lower triangle of the matrix has been omitted.
The two-electron energy expression parameters differ from the default settings, and are defined as:
x22 = 2p2 : x23 = 0 : x33 = 2q2 y22 = -p2 : y23 = pq : y33 = -q2
The general stategy within a given iteration cycle is to construct revised (hopefully improved molecular orbitals) as linear combinations of the input molecular orbitals to that cycle. Let w denote the row vector of trial molecular orbitals, and let w' denote a row vector of iterated molecular orbitals. Then:
w' = w(I + D) 7.4
where I denotes the identity matrix of order NEWBAS, and D denotes a NEWBAS*NEWBAS skew matrix, containing the molecular orbital mixing coefficients. The revised molecular orbitals are not orthonormal beyond first order. However, they are eventually rendered completely orthonormal by a subsequent application of the symmetric orthonormalisation scheme. Clearly what is required is a scheme for the computation of the matrix elements of D. For the mixing of two orbitals within a given shell, the corresponding matrix element of D is set to zero, since such mixing does not affect the energy. For mixing of two molecular orbitals (m and n) within different shells (one of which may be the virtual shell), the following formula is used:
Dmn = lMN (gMN / (-hmn - aMN)) 7.5 where gmn = (dE/dDmn)D = 0 7.6 hmn = (d2E/d2Dmn)D = 0 7.7
lMN : A user supplied input parameter, known as the damp factor. The same damp factor is used for all mixings between orbitals belonging to a given pair of shells. The default value for all damp factors is unity.
aMN : A user supplied input parameter, known as the level shifter. The same level shifter is used for all mixings between orbitals belonging to a given pair of shells. The default value for all level shifters is zero.
The lMN and aMN parameters may be user defined by means of the SHELLS directive. Before outlining the formulae used to evaluate the first and second derivatives of the energy, the Coulomb and exchange operators for a given shell via their matrix elements, are defined as follows:
JabM = [ / wm(2)wm(2) 1/r12 a(1)b(1) dt1 dt2 7.8
m
M
KabM = [ / wm(2)a(2) 1/r12 wm(1)b(1) dt1 dt2 7.9
m
M
further define operators as:
LM = WM F + 2 [ xMN JN + 2 [ yMN KN 7.10
N N
The type of molecular orbital mixing are considering is given by:
wm = wm + dnm wn + 1/2 dnm2 wm 7.11 wn = wn - dnm wm - 1/2 dnm2 wn 7.12
where one included the second order changes which will be induced by the symmetric orthonormalisation scheme. The following formulae define the first and second energy derivatives:
( dE )
------ = 2(LnmM - LnmN) 7.13
(dDnm)
( dE2 )
------- = 2(LnnM - LnnN - LmmM + LmmN
(dD2mn)
+ Jmn(2yMM + 2yNN - 4yMN)
+ Kmn(4xMM + 4xNN - 8xMN - 2yMM + 2yNN - 4yMN)) 7.14
To avoid the expense of computation of the Coulomb and exchange integrals, Jmn and Kmn which appear in equation 7.14, these are approximated as averages, thus:
Jmn = ( [ [ Jmn ) / (NM * NN) 7.15
m
M n
N
Kmn = ( [ [ Kmn ) / (NM * NN) 7.16
m
M n
N
where NM and NN denote the number of molecular orbitals in shells M and N respectively.
When calculating the D-matrix elements between orbitals of shells M and N, the program canonicalises these molecular orbitals over the operator (LM - LN), computes the appropriate D-matrix elements in this canonicalised basis, then transforms back to the original molecular orbital basis.
Since the energy is invariant to unitary mixing of molecular orbitals within a given shell, it is necessary to specify canonicalisations conditions before a unique definition of each molecular orbital is obtained. The procedure used by the program is to require, for a given shell (M), that a user defined Fock operator HM, is rendered diagonal in the basis of the molecular orbitals comprising the given shell. The general form of the Fock operator is given by:
HM = F + [ cMN JN + [ dMN KN 7.17
N N
where cMN and dMN are parameters which may be user specified by means of the JC and KC lines of the SHELLS directive, or are set to the default values cMN = WN, dMN = -WN /2.
Two molecular orbitals wm1 and wm2, within a given shell M, will therefore obey the following equations:
(wm1 !HM! wm1) = em1 7.18 (wm2 !HM! wm2) = em2 7.19 (wm1 !HM! wm2) = 0 (m1=/m2) 7.20
where em1 and em2 are referred to as the orbital energies of molecular orbitals wm1 and wm2 respectively. The user should note that the canonicalisations which the user specifies will not affect the rate of convergence of the iterative procedure. The procedure for determining canonicalisation has been suppiled simply so that the user may obtain molecular orbitals which have been rendered canonical in the most convenient form for the user own purposes.
The following examples give details of the energy expression parameters to be used with the GRHF program when studying linear molecule configurations of the form s1pn, n = 1 to 3. Such configurations may be coupled to a closed shell kernel of [+ symmetry, the closed shell having an occupation number (one-electron energy expression parameter) of 2. All two-electron energy expression parameters involving the closed shell (in interaction with itself or with open shells) take the default values as calculated by the GRHF program, and therefore need not be specified. In the following, two open shells are considered, the first being a single sigma orbital of occupation 1, the second consisting of a pair of degenerate pi orbitals occupied by n electrons, and having an occupation number of n/2.
Table 6.4 details the possible states derivable from the s1pn configuration. The detailed form of the wavefunctions for those components of highest ms values is also given.
Table 6.4
_________
Wavefunctions of the States Derived from s1pn
Configuration State Wavefunction
------------------------------------------------------------------------------
s1p1 1PI (!sx! - !sx!) / (2)**0.5
and (!sy! - !sy!) / (2)**0.5
3PI !sx!
and !sy!
s1p2 2[+ (!sxx! + !syy!) / (2)**0.5
2[- (2!sxy! - !sxy! + !syx!) / (6)**0.5
2D (!sxx! - !syy!) / (2)**0.5
and (!sxy! - !syx!) / (2)**0.5
4[- !sxy!
s1p3 1PI (!sxxy! - !sxxy!) / (2)**0.5
and (!syyx! - !syyx!) / (2)**0.5
3PI !sxxy!
and !syyx!
Table 6.5 details the Coulomb and exchange energy parameters for all the states derivable from the s1pn configurations. Note that the 4[- state derived from the s1p2 is of the half closed shell type, and may therefore also be treated using the RHF program, the open shell being s and both pi orbitals in that case. Where the values of the energy expression parameters are quoted in parentheses, they take the default value as calculated by the GRHF program, and need not therefore be specified through the SHELLS directive.
Table 6.5
_________
Coulomb and Exchange Energy Expression Parameters for the States
Derived from s1pn
! Coulomb ! exchange !
-----------------------------------------------------------------
Configuration ! State ! ss ! sp ! pp ! ss ! sp ! pp !
-----------------------------------------------------------------
s1p1 1PI (1/2) (1/4) 0 -1/2 1/4 0
3PI (1/2) (1/4) 0 -1/2 -1/4 0
s1p2 2[+ (1/2) (1/2) 0 -1/2 (-1/4) 0
2[- (1/2) (1/2) (1/2) -1/2 1/4 -1/2
2D (1/2) (1/2) 1/4 -1/2 (-1/4) 0
4[- (1/2) (1/2) (1/2) -1/2 -1/2 -1/2
s1p3 1PI (1/2) (3/4) 1 -1/2 0 -1/2
3PI (1/2) (3/4) 1 -1/2 -1/2 -1/2
Suppose that the first ten molecular orbitals comprise a closed shell kernal, orbital 11 is to be singly occupied s orbital, and orbitals 12 and 13 the triply occupied pi shell. The appropriate SHELLS directive for the 1PI state of the s1p3 configuration would look like:
SHELLS
1 TO 10
11
12 13
END
2 1 1.5
JE 3 3 1
KE 2 2 -0.5
KE 2 3 0
KE 3 3 -0.5
END
The following examples give details of energy expression parameters to be used with the GRHF program when studying linear molecule configurations of the form (p)1(p*)1, (p)3(p*)1 and (p)3(p*)3, where p and p* denote two orthogonal sets of p orbitals; these denote orbitals belonging to these sets by the symbols x, y, x* and y*. It is important, when using the energy expression parameters given in the following text, to ensure that the orbitals x and y (and similarly for x* and y*) remain properly parallel throughout the iterations. This gives rise to a difficulty: because the x and y orbitals are degenerate, so there is nothing in principle to prevent the GRHF program (or any other of the Hartree-Fock programs) from constructing the p orbitals as linear combinations of the x and y orbitals, and p* orbitals as non-parallel linear combinations of the x* and y* orbitals. Use of the ALIGN directive, allows the use to specify the alignment of such degenerate sets of orbitals.
The open shell systems considered may be coupled to a closed shell kernel of [+ symmetry, the closed shell will have a one-electron energy expression parameter of 2, and all two-electron energy expression parameters involving the closed shell (either in interaction with itself or with the open shells) take the default values as calculated by the GRHF program, and will not be specified here. The one-electron energy expression parameters for the p shells are either 1/2 or 3/2 depending upon whether they are singly or triply occupied respectively.
Table 6.6 details the possible states derivable from the (p)1(p*)1, (p)3(p*)1 and (p)3(p*)3. The form of the wavefunctionsfor those components of highest ms value is also given. The symbol P is used to denote the first part of the wavefunction (as written down) after permutation of the symbols x and y.
Table 6.6
_________
Configuration State Wavefunction
------------------------------------------------------------------------------
(p)1(p*)1 1[+ (!xx*! - !xx*! + P) / 2
1[- (!xy*! - !xy*! - P) / 2
1D (!xy*! - !xy*! + P) / 2
and (!xx*! - !xx*! - P) / 2
3[+ (!xx*! + P) / (2)**0.5
3[- (!xy*! - P) / (2)**0.5
3D (!xy*! + P) / (2)**0.5
and (!xx*! - P) / (2)**0.5
(p)3(p*)1 1[+ (!xxyy*! - !xxyy*! + P) / 2
1[- (!xxyx*! - !xxyx*! - P) / 2
1D (!xxyx*! - !xxyx*! + P) / 2
and (!xxyy*! - !xxyy*! - P) / 2
3[+ (!xxyy*! + P) / (2)**0.5
3[- (!xxyx*! - P) / (2)**0.5
3D (!xxyx*! + P) / (2)**0.5
and (!xxyy*! - P) / (2)**0.5
(p)3(p*)3 1[+ (!xxyx*x*y*! - !xxyx*x*y*! + P) / 2
1[- (!xxyy*y*x*! - !xxyy*y*x*! - P) / 2
1D (!xxyy*y*x*! - !xxyy*y*x*! + P) / 2
and (!xxyx*x*y*! - !xxyx*x*y*! - P) / 2
3[+ (!xxyx*x*y*! + P) / (2)**0.5
3[- (!xxyy*y*x*! - P) / (2)**0.5
3D (!xxyy*y*x*! + P) / (2)**0.5
and (!xxyx*x*y*! - P) / (2)**0.5
In the following it is assumed that 4 open shells have been defined, namely the x, y, x* and y* orbitals. The energy expression parameters for the xx, xy and yy interactions (and similarly for the x*x*, x*y* and y*y* interactions) depend only upon occupation number, and take the values 1 and -1/2 for Coulombic and exchange parameters respectively in the case of triply occupied p orbitals, and 0 for both Coulombic and exchange parameters in the case of singly occupied p orbitals. The parameters for the xx*, xy*, yy* and yx* interactions depend upon the particular state being studied, and are listed in table 6.7. Note that only the xx* and xy* parameters are listed explicitly, because yy* = xx*, yx* = xy*. In the case of the 3D state derived from (p)3(p*)1 it is possible to perform the calculation using only two open shells, namely the p and p* shells, consisting of the x, y, x* and y* orbitals respectively. In this case the Coulomb energy expression parameters of 1, 3/8 and 0 for the pp, pp* and p*p* interactions respectively, and -1/2, -1/4 and 0 for the pp, pp* and p*p* exchange parameters respectively. Note that, as in table 6.7, values quoted in parentheses take the default values as calculated by the GRHF program, and need not be specified in the SHELLS directive. Attention is drawn to the degeneracy of the 3[- and 1[- states of the (p)3(p*)1 configuration in the Hartree-Fock approximation. There is no need to do separate SCF calculations on these states. The results for the (p)3(p*)1 and (p)3(p*)3 configurations have been checked by comparison with the results of J.B. Rose and V. Mckoy [9], after appropriate conversion of these authors results to ATMOL standards. In all cases complete agreement was found.
Table 6.7
_________
Coulombic and Exchange Energy Expression Parameters for the States
Derived from (p)1(p*)1, (p)3(p*)1 and (p)3(p*)3 Configurations
! Coulombic Exchange !
----------------------------------------------------------
Configuration ! State ! xx* ! xy* ! xx* ! xy* !
----------------------------------------------------------
(p)1(p*)1 3[+ 3/8 -1/8 -1/4 -1/4
1[+ 3/8 -1/8 1/4 1/4
3[- -1/8 3/8 1/4 -3/4
1[- -1/8 3/8 -1/4 3/4
3D (1/8) (1/8) -1/4 1/4
1D (1/8) (1/8) -1/4 3/4
(p)3(p*)1 3[+ 1/8 5/8 -1/4 -1/4
1[+ 1/8 5/8 3/4 -5/4
3[- 5/8 1/8 -1/4 -1/4
1[- 5/8 1/8 -1/4 -1/4
3D (3/8) (3/8) -1/4 -1/4
1D (3/8) (3/8) -1/4 3/4
(p)3(p3*) 3[+ 11/8 7/8 -3/4 -3/4
1[+ 11/8 7/8 -1/4 -1/4
3[- 7/8 11/8 -1/4 -5/4
1[- 7/8 11/8 -3/4 1/4
3D (9/8) (9/8) -3/4 -1/4
1D (9/8) (9/8) -1/4 -3/4
Example 1:
Suppose that the molecular orbitals 9 and 10 comprise the x and y orbitals respectively, to be populated by 3 electrons, whilst orbitals 11 and 12 comprise the x* and y* orbitals, to be populated by 1 electron. Orbitals 1 to 8 comprise a closed shell [+ kernel. The appropriate SHELLS directive for the 1[+ state of the (p)3(p*)1 configuration would look like:
SHELLS
9
10
11
12
1 TO 8
END
1.5 1.5 0.5 2
JE 1 1 1
JE 1 2 1
JE 2 2 1 Defines the pp interaction
KE 1 1 -0.5
KE 1 2 -0.5
KE 2 2 -0.5
JE 3 3 0
JE 3 4 0
JE 4 4 0 Defines the p*p* interaction
KE 3 3 0
KE 3 4 0
KE 4 4 0
JE 1 3 0.125
JE 1 4 0.625
JE 2 3 0.625
JE 2 4 0.125 Defines the pp* interaction
KE 1 3 0.75
KE 1 4 -1.25
KE 2 3 -1.25
KE 2 4 0.75
Example 2
The appropriate SHELLS directive for the 3D state for the (p)3(p*)1 configuration, where one takes advantage of the possibility of using the only two shells, is given as:
SHELLS
9 10
11 12
1 TO 8
END
1.5 0.5 2
JE 1 1 1
KE 1 1 -0.5
JE 2 2 0
KE 2 2 0
JE 1 2 0.375
KE 1 2 -0.25
END
where the orbital structure is as given in Example 1. Note that the Coulombic energy expression parameter for the pp* interaction could be omitted, since it takes default values as calculated by the GRHF program.
The aim of this module (code named UHF), is the minimisation of the energy of spin unrestricted Hartree-Fock single determinental wavefunctions constructed using different spatial orbitals for different spins. The user should note that the NCORE and NOPEN parameters which appear on the first data line of data input are interpreted as follows:
NALPHA = NCORE+NOPEN
NBETA = NCORE
where NALPHA and NBETA denote the number of occupied alpha-spin and beta-spin molecular orbitals, respectively. The directives for the UHF module may be presented in any order, unless specifically requested.
See the TITLE directive of closed shell program.
See the MAXCYC directive of closed shell program.
See the MFILE directive of closed shell program.
See the SIZE directive of closed shell program.
See the CTRANS directive of closed shell program.
The general form of this directive is described in subparagraph 4.7. One extra parameter, namely ANIL, is allowed, and if used causes the program to output results concerning the effects of projection of the first contaminating spin multiplet from the UHF wavefunction.
See the PUNCH directive of the closed shell program. Note that the two sets of eigen vectors, population and eigen values will be 'punched'. The first set corresponds to the a-spin molecular orbitals, the second to the b-spin molecular orbitals.
See the FORMAT directive of closed shell program.
See the OUTPUT directive of closed shell program.
See the IPRINT directive of closed shell program.
See the ACCURACY directive of closed shell program.
See the TIME directive of closed shell program.
See subparagraph 4.14. Note that the occupied-virtual orbitals switching may occur within either the a or b-spin sets of molecular orbitals.
This directive consists of one data line, read to variables TEXT,DA1, EA1,DB1,EB1,IBRK,DA2,EA2,DB2,EB2 using format (A,4F,I,4F).
TEXT should be set to the character string DAMP.
DA1,EA1,DB1,EB1 the a-spin damp factor and level shifter are given by DA1 and EA1, respectively. While DB1 and EB1 are the b-spin damp and level shifters.
IBRK this designates the iterative cycle to which DA1,EA1,DB1 and EB1 hold too (inclusive of cycle IBRK).
DA2,EA2DB2,EB2 after cycle IBRK, these parameters will be used to specify the damp factors and level shifters.
An alternative form of the DAMP directive is permitted, consisting of only 5 data fields( ), read to variables TEXT,DA1,EA1,DB1 and EB1 using format (A,4F). This alternative form causes IBRK to be set to 999, DA2 and DB2 to be set to unity and EA2 and EB2 to zero. TEXT,DA1,EA1,DB1 and EB1 have their usual meanings.
The DAMP directive may be omitted, when the following default settings will be chosen:
DA1,DB1,DA2,DB2 = 1.0
EA1,EB1,EA2,EB2 = 0.0
IBRK = 999
The default settings correspond to performing the iterations in the manner prescribed by in reference [6].
See the LOCK directive of the open shell program. Note that if the LOCK directive is used, the program will issue an internal AUTO directive of the form:
AUTO 0
thus completely surpressing the automatic orbital switching feature.
Seven directives are available for obtaining trial eigen vectors as input to the first iterative cycle, namely VECTORS,START,ALPHAS,RESTORE, GETQ,EXTRA and COMBINE. The generic term LOADQ has been coined for use when referencing any of these directives. A LOADQ directive must appear in the input data, and each possibility is discussed in the following paragraphs (note that the syntax of the UHF-LOADQ directives is similar to the SCF-LOADQ directives).
The syntax of the first part of this directive is as for the corresponding closed shell module directive. This first stage of the input defines the trial alpha-spin eigen values. The second phase of the input is again a repeat of the closed shell module VECTORS syntax, and is used to define the trial beta-spin eigen vectors. The second phase of the directive may be omitted, when the trial beta-spin eigen vectors will be set equal to the trial alpha-spin eigen vectors.
See the START directive of the closed shell program. Notice that the trial a and b-spin eigenvectors will be equivalent if this directive is used.
See the ALPHAS directive of the closed shell program. The trial a and b-spin eigen vectors will be equivalent if this directive is used.
This directive consists of a single data line, read to variables TEXT, ISECTA,ISECTB using format (A,2I).
TEXT should be set to the character string RESTORE.
ISECTA should be set to identify the section on the DUMP FILE where the trial a-spin eigen vectors may be recovered. Such vectors will have normally been placed on the DUMP FILE by a previous run, of one of the other modules, the section used to store these eigen vectors being defined by means of the ENTER directive.
ISECTB should be used to specify the section of the DUMP FILE where the trial b-spin eigen vectors may be found. If ISECTB is omitted, then the trial b-spin vectors will be set equal to the a-spin vectors.
This directive causes restoration of trial eigen vectors from a foreign DUMP FILE, and consists of one data line read to variables TEXT,DDA, IBLKA,ISECTA,DDB,IBLKA,ISECTB using format (2A,2I,A,2I).
TEXT should be set to the character string GETQ.
DDA,IBLKA,ISECTA the trial a-spin vectors are to be found in the section specified by ISECTA of a foreign DUMP FILE residing on an ATMOL file assigned by the parameter DDA, starting at block IBLKA.
DDB,IBLKB,ISECTB the trial b-spin vectors are to be found in the section specified foreign DUMP FILE residing on an ATMOL section specified by ISECTB of a foreign DUMP FILE residing on an ATMOL file assigned by the DDB parameter, starting at block IBLKB. The DDB,IBLKB and ISECTB parameters may be omitted, when the trial b-spin vectors will be set equal to the trial a-spin vectors.
Example 1:
A UHF calculation has been performed at a geometry slightly different from the present case, and the DUMP FILE for this previous case is located to a data set assigned as ED4 starting at block 251. Assume that the converged a and b-spin vectors for this previous case are stored in sections 8 and 9 respectively of this DUMP FILE. The previous case vectors may be made available to the present run by means of the GETQ directive of the form:
GETQ ED4 251 8 ED4 9 251 9
If that a set of LCBF differs from that used in a previous calculation by the addition of a number of LCBF. Further, that these additional LCBF will not be of prime importance in the an updated calculation (for example, addition of polarisation functions).
The occupied a and b-spin orbitals of the previous calculation will provide a good starting point for the updated calculation, and the EXTRA directive has been implemented to be allow them to be used in this manner. The first data line is read to variables TEXT,DDA,IBLKA, ISECTA,DDB,IBLKB,ISECTB using format (2A,2I,A,2I).
TEXT should be set to the character string EXTRA.
DDA,IBLKA,ISECTA the DUMP FILE whereon the trial a-spin vectors (smaller basis) are to be found, resides on an ATMOL file assigned by the parameter DDA, and starting at the block specified by IBLKA. The trial a-spin vectors are to be found in the section specified by ISECTA.
DDB,IBLKB,ISECTB the DUMP FILE whereon the trial b-spin (smaller basis) vectors are to be found, resides on an ATMOL file assigned by the DDB parameter, starting at block IBLKB. The trial b-spin vectors are to be found in the section ISECTB.
Notice that the 'size' of the a- and b-spin trial eigen vectors must be the same (must be constructed from the same number of basis functions). DDB,IBLKB and ISECTB may be omitted, when the trial b-spin vectors will be set equivalent to the trial a-spin vectors.
Let the number of LCBF in the reference calculation be denoted by OLDBAS. Then the number of additional LCBF in the present calculation is given by NEWBAS-OLDBAS=NEXTRA. Subsequent data lines of the EXTRA directive will consist of NEXTRA integers read in free I-format, which are used to specify the LCBF which appear in the reference calculation. The last data field presented should contain the character string END. Note that it is assumed that those basis functions common to both present and reference calculations occur in the same order in both calculations.
The directive will issue an ORTHOG directive in SCHMIDT mode, so that trial eigen vectors will be automatically orthogonalised.
Example 1:
EXTRA ED4 1 1 ED4 1 2
11 12 13 END
specifies that trial a- and b-spin vectors are to be read from sections 1 and 2 respectively of the DUMP FILE stored in the data set assigned to ED4 and starting at block 1. LCBF 11 to 13 appear in the present calculation, yet did not appear in the reference calculation.
Example 2:
The use of the character string TO to abbreviated the method of specifying the list of extra LCBF is allowed. Thus the sequence:
EXTRA ED4 1 1 ED4 1 2
11 TO 13 END
is equivalent to that presented in example 1, above.
The syntax of the first part of this directive is as for the corresponding SCF module. This first stage of the input defines the a-spin vectors. The second phase of the input consists of a repeat of the SCF module COMBINE directive, and is used to define the trial b-spin vectors. This second phase of the directive may be omitted, when the trial beta vectors will be set equal to the trial a-spin vectors.
The notes on the LOADQ directive of the SCF module are equally applicable to the UHF module LOADQ directives.
See the SDIAG directive of closed shell program.
See the ORTHOG directive of closed shell program. Use of this directive causes orthonormalisation of trial alpha and beta-spin vectors.
See the NOORTH directive of closed shell program. Use of this directive causes no-orthonormalisation of trial alpha and beta-spin vectors.
The syntax of this directive is similar to that of the SCF module, except that the directive initiator is read to variables TEXT,SPIN using format (2A).
TEXT should be set to character string SWAP.
SPIN should be set to one of the character strings ALPHA or BETA, and will cause swapping of either alpha or beta spin vectors respectively. The SPIN parameter may be omitted (thus the directive is the same as the SCF module), and will cause a-spin vectors to be interchanged. Note that the first NALPHA and NBETA columns of the trial a- and b-spin vectors respectively will be deemed to be occupied.
This directive consists of a single data line, read to variables TEXT, ISECTA,ISECTB,ISECTC,ISECTD,ISECTE,ISECTF using format (A,6I).
TEXT should be set to the character string ENTER.
ISECTA specifies the section number of the DUMP FILE where the iterated a-spin vectors are to be placed.
ISECTB specifies the section number of the DUMP FILE where the iterated b-spin vectors are to be placed.
ISECTC specifies the section number of the DUMP FILE where the natural orbitals of the UHF wavefunction are to be placed. ISECTC may be set to zero, in which case the UHF natural orbitals will not be stored on the DUMP FILE.
ISECTD specifies the section number of the DUMP FILE where the spin natural orbitals of the UHF wavefunction are to be placed. ISECTD may be set to zero, in which case the UHF spin natural orbitals will not be on the DUMP FILE.
ISECTE specifies the section number of the DUMP FILE where the natural orbitals of the wavefunction produced by annihilation of the first contaminating spin multiplet from the UHF wavefunction are to be placed. ISECTE may be set to zero, in which case the natural orbitals of the annhilated UHF wavefunction are not produced if ANIL parameter of the FPRINT directive is not presented. Note also that the natural orbitals of the UHF function are identical before and after annihilation. The only difference lies in the occupation numbers.
ISECTF specifies the section number of the DUMP FILE where the spin natural orbitals of the annihilated UHF wavefunction are to be placed. ISECTF may be set to zero, in which case the spin natural orbitals of annihilated UHF function will not be routed to the DUMP FILE. Note that the spin natural orbitals of the annihilated UHF function are not produced if the ANIL parameter of the FPRINT directive is omitted. Note also that the spin natural orbitals of the UHF function are identical before and after the annihilation. The only difference lies in the occupation numbers.
The ENTER directive must be presented last, any directives presented after the ENTER directive will be ignored. The ENTER directive is illegal if a LOADQ directive has not been found in the preceding data.
This directive may be used as an alternative to the ENTER directive, and consists of a single data line, read to variables TEXT,ISECTA, ISECTB using format (A,2I).
TEXT should be set to the character string STOP.
ISECTA,ISECTB specifies the section numbers of the DUMP FILE where the trial a- and b- spin vectors are to be written.
The effect of the STOP directive is as described for the corresponding SCF module STOP directive.
The purpose of the UHF program is to minimise the energy of a single determinental wavefunction constructed from NALPHA orthonormal a-spin molecular orbitals, and NBETA orthonormal b-spin molecular orbitals. From NEWBAS linearly independent LCBF one may construct NVIRTA = NEWBAS - NALPHA virtual a-spin molecular orbitals, such that the complete set of a-spin molecular orbitals are orthonormal. Similarly, NVIRTB = NEWBAS - NBETA virtual b-spin molecular orbitals may be constructed. Let 0, X1a, X2a, X1b and X2b denote row vectors of the LCBF, occupied and virtual a-spin molecular orbitals, and occupied and virtual b-spin molecular orbitals respectively. Then:
(X1a:X2a) = 0(Q1a:Q2a) = 0Qa 8.1 (X1b:X2b) = 0(Q1b:Q2b) = 0Qb 8.2
where each column of Qa and Qb contains the coefficients of the LCBF within a given a- or b-spin molecular orbital respectively. Now introduce the a- and b-spin density matrices, Ra and Rb respectively:
Ra = Q1aQ1a+ 8.3 Rb = Q1bQ1b+ 8.4
Ra and Rb are invariant to unitary mixing between the a- and b-spin molecular orbitals. The total wavefunction is invariant to such mixing, being uniquely defined by Ra and Rb. The electronic energy expression is given by:
E = trRaF +trRbF + 1/2trRaJ[Ra] + 1/2trRbJ[Rb] + trRaJ[Rb]
- 1/2trRaK[Ra] - 1/2trRbK[Rb] 8.5
where F is the matrix representative of the one-electron Hamiltonian operator, J[R] and K[R] denote Coulomb and exchange matrices respectively, as described in subparagraph 4.xx. Now define:
R = Ra + Rb 8.6 P = (Ra - Rb) / 2 8.7
where R is the total density matrix of the system, while P is the spin
G = J[R] - 1/2K[R] 8.8 Z = K[P] 8.9 H = F + G 8.10
The energy expression may be expressed in the computationally convenient form:
E = 1/2trR(F+H) - trPZ 8.11
It is convenient to define the a- and b-spin Fock matrices, Ha and Hb:
Ha = H - Z 8.13 Hb = H + Z 8.14
Physically significant variations in the molecular orbitals as follows:
(a) Mixing of a-spin occupied with virtual molecular orbital. Small such variations may be written:
w1a = (w1a:w2a) !Ia!
!--! 8.14
!Da!
where Ia is an identity matrix of the order NALPHA, and Da is a NVIRTA*NALPHA matrix of mixing coefficients.
(b) Mixing of b-spin occupied with virtual molecular orbitals. Small such variations may be written:
w1b = (w1b:w2b) !Ib!
!--! 8.15
!Db!
where Ib is an identity matrix of order NBETA, and Db is a NVIRTB*NBETA matrix of mixing coefficients.
The perturbed molecular orbitals are not orthonormal beyond first order. Fortunately, the energy expression is required to be correct to first order only, and is given by:
NVIRTA NALPHA NVIRTB NBETA
E = E0 + 2 [ [ Dkia Hkia + 2 [ [ Dljb Hljb
k i l j
+ higher terms
where Hkia denotes the matrix element connecting the k'th virtual with the i'th occupied a-spin orbital over the a-spin Fock operator, Ha, and Hljb denotes the matrix element connecting the l'th virtual with the j'th occupied b-spin orbital over the b-spin Fock operator, Hb. It is clear that the conditions for self consistency are satisfied when all Fock matrix elements, Hkia and Hljb are zero. Such that:
( dE ) ------ = 2Hkia 8.17 (dDkia)D=0 ( dE ) ------- = 2Hljb 8.18 (dDljb)D=0
The form of equation 8.16 indicates that any minimisation procedure capable of producing values for the elements of the D-matrices such that all Dkia are of opposite sign to the corresponding Hkia, and all Dljb are of opposite sign to the corresponding Hljb, and such that the D-matrix elements can be chosen sufficiently small in magnitude that the higher terms of equation 8.16 are smaller in magnitude than the first order terms, will be capable of generating an iterated wavefunction of lower energy than the trial wavefunction, and therefore will guarantee convergence. Note that convergence to the lowest energy stationary point will not necessarily be guaranteed. The ATMOL procedure is as follows:
First construct the Fock operators Ha and Hb in the basis of the trial a- and b-spin molecular orbitals, then incorporate damp factors and level shifters to produce modified Fock operators, HMODa and HMODb as follows:
a !(Ha)oo :la(Ha)ov !
H = !--------------------! 8.19
MOD !la(Ha)vo:(Ha)vv+ubI1!
and
b !(Hb)oo :la(Hb)ov !
H = !--------------------! 8.20
MOD !lb(Hb)vo:(Hb)vv+ubI2!
where:
(Ha)oo and (Hb)oo: denote a block of Ha and Hb in the basis of the occupied a- and b-spin trial molecular orbitals.
(Ha)vv and (Hb)vv: denote a block of Ha and Hb in the basis of the virtual a- and b-spin molecular orbitals.
(Ha)ov and (Hb)ov: denote a block of Ha and Hb connecting all a-spin occupied with virtual molecular orbitals and all b-spin occupied with virtual molecular orbitals.
ua and ub: denote a- and b-spin level shifters.
la and lb: denote a- and b-spin damp factors.
Both Fock operators are diagonalised, and unless a LOCK directive is used, the resulting eigen vectors are ordered according to the aufbau principle based on the associated level shifted eigen values. If the LOCK directive is used, ordering of the eigen value columns is based on the principle of maximum overlap of the iterated vector with the trial vector. Whichever ordering scheme is adopted, the first NALPHA columns of the a-spin vector array define iterated occupied a-spin occupied molecular orbitals as linear combinations of the trial a-spin molecular orbitals. Similarly, the first NBETA columns of the b-spin vector array define iterated occupied b-spin molecular orbitals.
Analysis of the above diagonalisation procedure shows that if the damp parameters are set to the default value of unity, and the LOCK directive is not used, convergence in a downward direction on the energy surface can be guaranteed if a sufficiently large and positive valued level shifter is employed. The analysis is closely similar to that outlined for the closed shell case. Note that the choice of unity for the damp factor and zero for the level shifters corresponds to performing the iterations in the manner prescribed by J.A. Pople and R.K. Nesbet [10].
In the case that NALPHA=NBETA, it is easy to show that a wavefunction whose energy is stationary on the closed shell RHF surface will also have a stationary energy on the UHF surface. Furthermore, it is usually found that for closed shell molecules near their equilibrium geometry the minimum energy solutions of the RHF and UHF equations are identical. However if bond lengths much greater than the equilibrium are employed, then the RHF solution corresponds only to a saddle point on the UHF surface, so that there exists a lower energy stationary point on the UHF surface. Suppose a trial wavefunction for a UHF closed shell case, has identical a- and b-spin molecular orbitals. Unless special precautions are taken, the present method, which treats a- and b-spin molecular orbitals symmetrically, will converge to a self consistent UHF wavefunction which is identical to the corresponding RHF wavefunction. In such cases it is wise for the user to force an asymmetry of the treatment of the a- and b-spin orbitals. Such asymmetry may be achieved in two ways:
(a) The use of different values for the a- and b-spin damp factors (la and lb) or the level shifters (ua and ub), will produce asymmetry.
(b) Dissimilar trial a- and b-spin vectors may be generated by performing a UHF calculation on some closely related open shell species. The vectors generated may be made available to the closed shell UHF case by means of a RESTORE directive. Since the purpose of the open shell calculation is simply to provide dissimilar a- and b-spin vectors, there is no need to procede to convergence, only one cycle would be sufficient.
As is known, the UHF wavefunction is not in general an eigen function of the S2 operator. The wavefunction can be spectrally analysed into components which are eigen functions of the S2 operator, such components having spin quantum numbers, s, which lie in the range:
1/2(NALPHA + NBETA) > s > 1/2(NALPHA - NBETA)
Each component is a of multiplicity 2s+1, with eigen values (with respect to S2) of s(s+1). Assume that the UHF wavefunction consists predominantly of the component of lowest possible spin quantum number, m, where:
m = 1/2(NAPLHA - NBETA)
and that the principal contamiant is that multiplet of spin quantum number m+1. The operator:
As = S2 - s(s+1)
when acting upon the UHF wavefunction will annihilate the component of spin quantum number s. Clearly, the operator Am+1 will annihilate the contaminant of spin quantum number m+1. The wavefunction:
WAUHF = N Am+1 wUHF
will be referred as the annihilated UHF wavefunction (N is a normalising function). In the implementation of the above annihilation procedure, the form of the effective S2 operator takes, is:
S2 = [' Pij + 1/4[(na - nb)2 + 2na + 2nb]
ij
where Pij interchanges spin functions of electrons i and j, the prime indicates on the summation that only interchanges involving differing spin functions are to be summed. Hence:
Am+1 = [' Pij + 2NBETA - NALPHA - 2
ij
= [' Pij + NBETA - 2(m+1)
ij
The theory necessary for the analysis of WAUHF can be found in paper by T. Amos and L.C. Snyder [11].
The purpose of this module (code named BOYS), is the localisation of molecular orbitals according to the criterion of Foster and Boys [12].
The user should note that NCORE and NOPEN parameters which appear on the standard first data line have no significance for the BOYS module, although they are submitted to certain standard checking procedures, and its is advised that both these parameters should be set to unity. The directives of the BOYS module following the first data line can be presented in any order, unless previously stated.
See the TITLE directive of closed shell program.
See the MAXCYC directive of closed shell program.
See the SIZE directive of closed shell program.
See the CTRANS directive of closed shell program.
See the FPRINT directive of closed shell program. Note that only the parameters S,T,T+V,X,Y and Z are relevant.
See the PUNCH directive of closed shell program. Note that only the parameter VECTORS is relevant.
See the FORMAT directive of the closed shell program. Note that only the FORM1 parameter is relevant.
See the OUTPUT directive of closed shell program.
This directive consists of one data line with the character string IPRINT in the first data field. If used, printing of the dipole centroids of the molecular orbitals involved in the localisation process will be produced from each iterative cycle.
See the ACCURACY directive of closed shell program.
This directive is used to define those molecular orbitals which take part in the localisation process. The first data field consists of the character string ACTIVE in the first data field. Subsequent data lines are read to an array (IACTIV(I),I=1,NACT) using free I-format. The last data field presented should be the character string END.
Example 1:
ACTIVE
2 3 4 5 6 7 9 END
Declares molecular orbitals 2 to 7 inclusive, plus molecular orbital 9 to be active in the localisation process.
Example 2:
The sequence
ACTIVE
2 3 4 5 6 7
9
END
has an equivalent effect to that of example 1.
Example 3:
An abbreviated form of specifying the list of active molecular orbitals is to invoke the character string TO, to link together a sequence of consecutive numbered active molecular orbitals.
ACTIVE
2 TO 7 9 END
The above sequence is equivalent to examples 1 and 2.
Example 4:
ACTIVE
2 TO 7 9 TO 14
END
specifies that molecular orbitals 2 to 7, and 9 to 14 are to be active in the localisation process.
This directive consists of one data line, with the character string LOCK in the first data field. In the presence of the LOCK directive, the program will seek a stationary value for the functional, and proceed in the direction of 'minimum change'. The use of the LOCK directive may cause the procedure to fail to converge. In the absence of the LOCK directive, the program will unconditionally maximise the functional defined in equation 5.
Seven directives are available for obtaining the molecular orbitals to be localised; namely VECTORS,START,ALPHAS,RESTORE,GETQ,EXTRA and COMBINE. One LOADQ directive must appear in the input data, see SCF LOADQ directives for syntax.
It is expected that the RESTORE directive will be the most commonly used, either to retrieve molecular orbitals produced by one of the Hartree-Fock modules, or to retrieve partially localised orbitals produced by a run of the BOYS localisation module which was given insufficient iterations to achieve convergence.
See the SDIAG directive of closed shell program.
See the ORTHOG directive of closed shell program.
See the NOORTH directive of closed shell program.
See the SWAP directive of closed shell program. Note that the form of the ACTIVE directive requires the user to know in which order the molecular orbitals are to be found, so it is not expected that the SWAP directive would be heavily used in the BOYS module.
See the ORBSYM directive of closed shell program. Again it is not expected the use of the ORBSYM directive will be extensively used in the BOYS module, for the same reason given in the SWAP directive.
See the ALIGN directive of closed shell program.
See the MOPOPS directive of closed shell program.
See the EVALUES directive of closed shell program.
See the ENTER directive of closed shell program.
See the STOP directive of closed shell program.
The method of Foster and Boys generates localised orbitals by specifying that unitary transformation of the orbitals which maximises the sum of the orbital self quadratic repulsion integrals:
NACT
I = [ 9.1
i
This procedure may be shown equivalent to one where one maximises the sum of the squares of the distances between the dipole centroids of the active orbitals.
NACT 2
J = [ ( - ) 9.2
i>j
Equivalently one may maximise the sum of the squares of the distances of the dipole centroids of the active orbitals, measured from an artibrary point, in the present case the origin as specified in the data input for the molecular integrals run:
NACT 2
K = [ 9.6
i
The maximisation of K (equation 9.6) is accomplished by an iterative procedure, involving successive rotations of pairs of active molecular orbitals. One cycle of the iterative procedure is deemed to be complete when all possible such pairs have been subjected to a rotation.
Consider the two dimensional case involving active orbitals wk and wl. Express the rotations by the usual form of equations:
wk' = cosa wk + sina wl wl' = -sina wk + cosa wl
The functional, K, may be written as a function of a:
K(a) = K(0) - A + Acos4a + Bsin4a 9.4
where
A = ( - )2 / 4 - 2 B = ( - )2 . 2
hence
dK(a)
----- = - Asin4a + Bcos4a
da
Require solutions for a, such that the above equation = 0. Two solutions can be found.
Solution 1
a is chosen so that the following equations are satisfied:
sin4a = B / (A2 + B2)**0.5 cos4a = A / (A2 + B2)**0.5
This solution corresponds to maximising the functional with respect to a, and will be chosen if the LOCK directive is not requested. The functional acquires the value:
K(amax) = K(0) - A + (A2 + B2)**0.5
This solution also corresponds to the minimum change stationary point (solution with the smallest value of a and sina), in the event that A is positive, this will be used if the LOCK directive is presented.
Solution 2
a is chosen so that the following equations are satisfied:
sin4a = -B / (A2 + B2)**0.5 cos4a = -A / (A2 + B2)**0.5
This solution corresponds to minimising the functional with respect to a, and would not be chosen if the LOCK directive was not presented. The functional acquires the value:
K(amin) = K(0) - A - (A2 + B2)**0.5
This solution corresponds to the minimum change solution, in the event of A being negative, and will be chosen if the the LOCK directive is presented.
It is obvious that if Solution 1 is used unconditionally in the case that the LOCK directive was not used, or choose between the two possible Solutions by monitoring the sign of A in the case where a LOCK directive was used. One is left with equations of the form:
sin4a = 'a' cos4a = 'b'
to solve (where a2 + b2 =1). The procedure is, determine:
c = sin-1 a
This yields a value for c such that -pi/2 < c < pi/2, given a suitable definition of the sin-1 function. If b is positive, a may be determined from:
a = c / 4
If b is negative, then a may be determined from either:
a = (pi - c) / 4
in the case that c is positive valued, or:
a = - (pi + c) / 4
in the case that c is negative valued.
The above procedure is most numerically stable when determining values for a which are small, this case being by far the most frequently met in the practical application.
The possible ATMOL error codes with a brief explanation are given in the following table:
Error Code Explanation
__________ ___________
14 Directive should not be used with the selected module.
15 Illegal value for NBASIS (1 < NABSIS < 256).
16 Directive unknown.
31 Illegal value for NCORE or NOPEN.
32 An attempt has been made to use an undefined MAIN FILE.
33 The MFILE directive to specify the MAIN FILE is missing.
34 The parameters specifying the symmetry of the LCBF are
not valid. May occur in the ALIGN directive.
35 Illegal value for ITRAN (in CTRANS directive).
36 Illegal value for NTRAN (in CTRANS directive).
37 Illegal value for NEWBAS (1 < NEWBAS < NBASIS+1).
38 The SWAP directive appears before a LOADQ directive.
39 The parameters of the SWAP directive are invalid.
40 An ENTER or STOP directive appears before a
LOADQ directive.
41 Two LCBF which have been given different symmetry tags
in the ALIGN directive are not orthogonal.
42 AFN in MFILE, GETQ, EXTRA or COMBINE
directive not known.
43 The ALIGN directive appears before a LOADQ directive.
46 The ORBSYM directive appears before a LOADQ directive.
47 The parameters specifying the symmetry of the LCBF
are not valid. May occur in the ORBSYM directive.
48 A molecular orbital of the specified symmetry was
not found during processing of the ORBSYM directive.
49 Invalid sequence of parameters specifying the required
order of molecular orbitals in the ORBSYM directive.
50 Invalid parameter in WIDTH pre-directive.
51 The MOPOPS directive appears before a LOADQ directive.
52 The EVALUES directive appears before a LOADQ directive.
61 Index block of DUMP FILE not in correct format.
62 ATMOL block with invalid checksum has been read,
or input/output error on ATMOL file. If the
latter, a finite VSOS error code will be given
whose explanation will be found in [8].
63 An attempt has been made to use a Section number
on the DUMP FILE outside the allowed range.
(1 < Section number < 190).
64 An attempt has been made to retrieve an undefined
Section from the DUMP FILE.
65 A retrieved Section on the DUMP FILE is of the
wrong TYPE.
66 ATMOL data set not assigned.
67 Illegal search of ATMOL data set.
68 A data field was read in F-format, and an illegal
character was detected.
69 A data field was read in I-format, and an illegal
character was detected.
70 The SIZE directive specified a maximum size less than
than the current size of the DUMP FILE.
71 An attempt has been made to expand the DUMP FILE
beyond its maximum size (as specified by a SIZE
directive).
72 An attempt has been made to overwrite a Section on
the DUMP FILE with a Section of greater length.
92 The absolute value of a coefficient in the CTRANS
directive is less than 10**(-8).
99 Module not known. This error may occur when
attempting to interpret the first parameter of the
first data line.
150 The eigen vectors retrieved by either of the RESTORE,
GETQ, EXTRA or COMBINE directives are invalid.
151 The CTRANS data in the input stream is incompatible
with that retrieved from the DUMP FILE during execution
of a RESTORE directive.
153 An illegal parameter has been detected in an EXTRA
directive.
154 Incorrect number of EXTRA LCBF specified.
156 LCBF declared twice in the EXTRA directive.
159 Invalid value for NELEC in the ENERGY directive.
160 Invalid parameters specified in the ACTIVE directive.
161 Invalid number of active molecular orbitals.
162 A molecular orbital has been specified more than once
in the ACTIVE directive.
163 Invalid molecular orbital declared to be a member of
a shell. May occur in the SHELLS directive.
164 Molecular orbital referenced more than once in a SHELLS
directive.
165 Invalid number of molecular orbitals in a shell.
166 Invalid number of shells.
167 Invalid reference to a shell on a JE, KE, JC, KC,
DAMP or SHIFT line of the SHELLS directive.
168 One of the character strings JE, KE, JC, KC, DAMP
or SHIFT was expected, but was not found. This
diagnostic may occur during processing of the
SHELLS directive.
171 An error has been detected in the data declaring
the origins of the LCBF of the AB system in the
COMBINE directive.
172 An error has been detected in the data declaring
the origins of the molecular orbitals of the AB
system in the COMBINE directive.
This job illustrates how to run the closed shell module of the SCF program. The Gaussian Integrals used in this example are generated via INTEGV [X]. The trial test MO's are generated via the START directive. The SCF vectors are placed in section 1 of the DUMP FILE, while the accuracy of the energy is set at 10**(-8) if convergence is achieved.
/*JOB JOBNAME,ACCOUNT,ST=(C20,LP=1,WS=256),PW=PASSWORD,TI=30,C=B
PATTACH,ATMOL.
ATTACH,ED2V,ED3V,ACC=RW.
REQUEST,ED7,RT=U.
SCF.
####S
LPAGE 1
CHANGE ED2 ED2V ED3 ED3V
SCF 25 5 0 1 ED3
TITLE
(H2O) SCF
MFILE
ED2
1
0
OUTPUT
DAMP 1 1 10 .5 .5
ACCURACY 1 8
MAXCYC 100
START
ENTER 1
####S
/*JOB JOBNAME,ACCOUNT,ST=(C20,LP=2,WS=512),PW=PASSWORD,TI=100,C=B
PATTACH,ATMOL.
ATTACH,DIMED2,DIMED3,ACC=RW.
REQUEST,ED7,RT=U.
SCF.
####S
LPAGE 2
FILE ED2 DIMED2 ED3 DIMED3
SCF 84 10 0 1 ED3
TITLE
(H2O)2 SCF
MFILE
ED2
1
0
OUTPUT
DAMP 1 2 8 .4 .4
MAXCYC 50
START
ENTER 1
####S
/*JOB JOBNAME,ACCOUNT,ST=(C20,LP=1,WS=256),PW=PASSWORD,TI=30,C=B
PATTACH,ATMOL.
ATTACH,ED2V,ED3V,ACC=RW.
REQUEST,ED7,RT=U.
SCF.
####S
LPAGE 1
FILE ED2 ED2V ED3 ED3V
RHF 25 4 1 1 ED3
TITLE
(H2O)+ RHF
MFILE
ED2
1
0
OUTPUT
D12 .9
D13 .89
D23 .91
MAXCYC 42
ACCURACY 3 5
LOCK
LEVEL .2 .2 6 .1 .1
RESTORE 1
ENTER 2
####S
/*JOB JOBNAME,ACCOUNT,ST=(C20,LP=1,WS=256),PW=PASSWORD,TI=30,C=B
PATTACH,ATMOL.
ATTACH,ED2V,ED3V,ACC=RW.
REQUEST,ED7,RT=U.
SCF.
####S
LPAGE 1
CHANGE ED2 ED2V ED3 ED3V
SERHF 25 3 2 1 ED3
TITLE
(H2O)++ SERHF
MFILE
ED2
1
0
OUTPUT
IPRINT
D12 .5
D13 .69
D23 .61
MAXCYC 20
ACCURACY 4 5
LOCK
ENERGY 2 0 -0.333333
LEVEL .2 .2 6 .1 .1
RESTORE 1
SWAP
3 4
END
ENTER 3
####S
/*JOB JOBNAME,ACCOUNT,ST=(C20,LP=1,WS=256),PW=PASSWORD,TI=30,C=B
PATTACH,ATMOL.
ATTACH,ED2V,ED3V,ACC=RW.
REQUEST,ED7,RT=U.
SCF.
####S
LPAGE 1
FILE ED2 ED2V ED3 ED3V
GRHF 25 5 1 1 ED3
TITLE
(H2O) GRHF
MFILE
ED2
1
0
OUTPUT
MAXCYC 50
SHELLS
1 TO 3
4
5
END
2 1.75 1.25
SHIFT 1 2 1
SHIFT 1 3 1
SHIFT 1 4 1
SHIFT 2 3 1
SHIFT 2 4 1
SHIFT 3 4 1
END
REST 1
ENTER 4
####S
/*JOB JOBNAME,ACCOUNT,ST=(C20,LP=1,WS=256),PW=PASSWORD,TI=30,C=B
PATTACH,ATMOL.
ATTACH,ED2V,ED3V,ACC=RW.
REQUEST,ED7,RT=U.
SCF.
####S
LPAGE 1
FILE ED2 ED2V ED3 ED3V
UHF 25 4 1 1 ED3
TITLE
(H2O)+ UHF
MFILE
ED2
1
0
OUTPUT
MAXCYC 42
ACCURACY 3 5
LOCK
LEVEL .2 .2 6 .1 .1
RESTORE 1
ENTER 5 6
####S
/*JOB JOBNAME,ACCOUNT,ST=(C20,LP=1,WS=256),PW=PASSWORD,TI=30,C=B
PATTACH,ATMOL.
ATTACH,ED2V,ED3V,ACC=RW.
REQUEST,ED7,RT=U.
SCF.
####S
CHANGE ED2 ED2V ED3 ED3V
LPAGE 1
BOYS 25 1 1 1 ED3
TITLE
(H2O) BOYS
MFILE
ED2
1
0
OUTPUT
ACTIVE
2 TO 5
END
MAXCYC 45
ACCU 1 5
RESTORE 1
ENTER 8
####S
[1] D. Moncrieff and V.R. Saunders, ATMOL-Introductory Notes.
[2] D.Moncrieff and V.R. Saunders, ATMOL-Integrals program.
[3] D. Moncrieff and V.R. Saunders, ATMOL-SERVICE program.
[4] Roothaan, C. C. Self-consistent field theory for open shells of electronic systems, J. Rev. Mod. Phys., 32, 179-185 (1960).
[5] V. R. Saunders and I. H. Hillier, A 'Level-Shifting' method for converging closed shell Hartree-Fock wave functions, Int. J. Quantum Chem. 7
[6] Jacobi diagonalisation ??
[7] R.F. Stewart, J. Chem. Phys.,431,52,(1970).
[8] M. F. Guest and V. R. Saunders, Mol. Phys. 28, 819 (1974). On methods for converging open-shell Hartree-Fock wave-functions
[9] J.B. Rose and V. Mckoy, J Chem Phys, 55, 1457 (1971)
[10] J.A. Pople and R.K. Nesbet, 'Self-Consistent Orbitals for Radicals', Journal of Chemical Physics, Vol. 22, No. 3, 1954, pp. 571-572
[11] T. Amos and L.C. Snyder, Unrestricted Hartree-Fock Calculations. I. An Improved Method of Computing Spin Properties, J. Chem. Phys. 41, 1773 (1964)
[12] Foster, J. M.; Boys, S. F. (1960). "Canonical Configuration Interaction Procedure". Reviews of Modern Physics. 32 (2): 300-302