13 Analysing the Wavefunction

GAMESS-UK includes a variety of tools for analysing wavefunctions, driven by the RUNTYPE ANALYSE directive. It is now possible to:

• Calculate a variety of 1-electron properties.
• Generate a localised orbital representation of an SCF wavefunction using either the dipole-centroid technique due to Foster and Boys, or the overlap-based criterion due to Pipek and Mezey [31].
• Provide graphical analysis of molecular wavefunctions. The program is capable of generating contour and perspective plots which depict:
  • the electron density associated with one or more molecular orbitals;
  • the amplitude of a molecular orbital;
  • a comparison of the density distribution in two or more molecular systems;
  • the interaction energy between a molecular distribution and a hypothetical point charge, generating the so-called electrostatic potential plot.
  Two types of plot may be generated by the program to provide a pictorial representation of a given density or potential function in a specified molecular plane:
  • a contour plot, with contours representing lines of constant value, depicting the spatial characteristics of the given function;
  • a perspective plot, with the values of the function in a given plane, displayed as a 3-D perspective picture.
  Both types of plot are generated from a grid of function values produced by the program.
• Perform a distributed multipole analysis (DMA) of an SCF wavefunction [32].
• Perform a more detailed Mulliken analysis, including both bond and orbital properties.

The user should note the following;

• While the various SCF modules provide default sections for eigenvector information, it will be necessary in the ANALYSE modules to specify, via the VECTORS directive, the specific eigenvectors to be analysed.
• Each mode of analysis typically requires one or more directives to specify the particular tasks required.

At present we restrict ourselves to sample data files for property evaluation, localised orbital analysis, graphical analysis, DMA and extended-Mulliken analysis. In each case we assume that the closed-shell SCF calculation on formaldehyde (I) has been successfully completed, and perform the requested analysis based on the SCF-MOs, as written to the Section 1 of the Dumpfile.

13.1 One-electron Property Evaluation

The following data sequence would be required in evaluating the electric field gradient at the carbon and oxygen nuclei.

          RESTART
          TITLE
          H2CO - 3-21G DEFAULT BASIS - 1-E PROPERTIES
          ZMATRIX ANGSTROM
          C
          O 1 1.203
          H 1 1.099 2 121.8
          H 1 1.099 2 121.8 3 180.0
          END
          RUNTYPE ANALYSE
          PROPERTY
          4 C
          4 O
          END
          VECTORS 1
          ENTER

Each one-electron operator is known to the user by an operator number; a full list of the available operators and associated numbers in given in Table 6 . The user specifies, under control of the PROPERTY directive, those properties to be be computed at any of the nuclei known to the system, by virtue of the TAGs defined in the z-matrix.

Table: The One-electron Operators and Operator Numbers

Operator	Operator	Operator	Operator
Number		Number
1	Potential	11	Third Moment (combined)
2	Diamagnetic Shielding	12	Hexadecapole Moment
3	Electric Field	13	Fourth Moment (even)
4	Electric Field Gradient	14	Fourth Moment (odd)
5	Dipole Moment	15	Overlap
6	Quadrupole Moment	16	Planar Charge Density
7	Diamagnetic Susceptibility	17	Line Charge Density
8	Second Moment	18	Charge Density
9	Octupole Moment	19	Isotropic ESR Coupling Constants
10	Third Moment	20	Anisotropic ESR Coupling Constants

The example above typifies the case where a single set of MOs are associated with the particular SCFTYPE, and as such may be input under control of the VECTORS directive to the properties package. A somewhat different approach is required when computing the one-electron properties derived from a wavefunction with more than one set of MOs (e.g., a UHF wavefunction), or in cases where only the total density matrix, and not an associated set of MOs, is available (e.g., in an MP2 calculation). In both cases, the user will need to generate the associated set of spinfree natural orbitals, and present these as input to the analysis module. Such orbitals are generated under control of the NATORB directive, which may used to route the natural orbitals to a nominated section on the Dumpfile.

The following data sequences would be required when evaluating the properties based on a UHF wavefunction. First, the data for the UHF calculation itself:

          TITLE
          H2CO - 3A2 UHF - 3-21G DEFAULT BASIS 
          MULT 3
          ZMATRIX ANGSTROM
          C
          O 1 1.203
          H 1 1.099 2 121.8
          H 1 1.099 2 121.8 3 180.0
          END
          SCFTYPE UHF
          NATORB 10 PRINT
          ENTER

Having routed the spinfree natural orbitals to section 10 on the Dumpfile, the properties calculation proceeds by nominating this section on the VECTORS line, thus:

          RESTART NEW
          TITLE
          H2CO - 3A2 UHF - 3-21G DEFAULT BASIS - 1-E PROPERTIES
          ZMATRIX ANGSTROM
          C
          O 1 1.203
          H 1 1.099 2 121.8
          H 1 1.099 2 121.8 3 180.0
          END
          RUNTYPE ANALYSE
          SCFTYPE UHF
          PROPERTY
          4 C
          4 O
          END
          VECTORS 10
          ENTER

As presented above, the NATORB directive will request generation of the spinfree natural orbitals. Two variants of the directive allow for (i) generation of the spin natural orbitals, and (ii) annihilation of the UHF wavefunction and subsequent generation of both spin- and spinfree NOs. The associated data requirements are straightforward:

• Generation of the spin NOs is driven by the keyword SPIN presented immediately after the NATORB initiator. Thus the data line

          NATORB SPIN 11 PRINT
  

  would result in the routing of the spin NOs to section 11 of the Dumpfile.
• Annihilation of the UHF wavefunction, and subsequent generation of the NOs from the annihilated density matrices is driven by specification of the keyword ANNIHILATE as the final character string on the NATORB data line. Thus the data sequence:

          NATORB 12 PRINT ANNIHILATE
  

  would route the spinfree NOs of the annihilated UHF wavefunction, AUHF, to section 12 of the Dumpfile.

The theory behind the AUHF analysis can be found in [33]. Note that the NOs of the UHF and AUHF wave function are in fact identical, the only difference lying in the occupation numbers.

Now let us consider the date requirements when computing properties at the optimum geometry derived from an MP2 calculation. First, the data for the MP2 geometry optimisation, where the spinfree natural orbitals at the optimised geometry are to be routed to section 20.

          TITLE
          H2CO - X1A1 - MP2 (DZ BASIS)
          ZMATRIX ANGSTROM
          C
          O 1 CO
          H 1 CH 2 HCO
          H 1 CH 2 HCO 3 180.0
          VARIABLES
          CO 1.203\CH 1.099\HCO 121.8
          END
          BASIS DZ
          RUNTYPE OPTIMISE
          SCFTYPE MP2
          NATORB 20 PRINT
          ENTER

Having routed the spinfree natural orbitals to section 20 on the Dumpfile, the properties calculation proceeds by nominating this section on the VECTORS line, thus:

          RESTART 
          TITLE
          H2CO - X1A1 - MP2/DZ BASIS - 1-E PROPERTIES
          ZMATRIX ANGSTROM
          C
          O 1 CO
          H 1 CH 2 HCO
          H 1 CH 2 HCO 3 180.0
          VARIABLES
          CO 1.203\CH 1.099\HCO 121.8
          END
          BASIS DZ
          RUNTYPE ANALYSE
          SCFTYPE MP2
          PROPERTY
          4 C
          4 O
          END
          VECTORS 20
          ENTER

Note the use of RESTART in restoring the optimized geometry from the Dumpfile.

13.2 Simplified Property Specification

In the examples above we have assumed that property evaluation is to be conducted under control of RUNTYPE ANALYSE, with explicit specification of the required one-electron properties. A simplified mechanism for property evaluation can be requested through presenting the data line

         PROPERTY ATOMS

after RUNTYPE and SCFTYPE specification. This will result in the default wavefunction analysis conducted after RUNTYPE processing being augmented with the computation of certain one-electron properties. The following points should be noted:

• the properties evaluated include the electrostatic potential, electric field, electric field gradient, and electron density at each of the atomic centres, plus the dipole, second moment, quadrupole moment, third and octupole moments, at the computed centre of mass of the system under study. In addition the spin densities will also be computed in the case of open shell systems.
• this analysis, if requested, is available on completion of SCF, OPTIMIZE, OPTXYZ, SADDLE, and CI processing.

The following data sequence would be required to generate the above list of properties on completion of an SCF calculation of the formaldehyde molecule.

          TITLE
          H2CO - 3-21G BASIS - SCF + DEFAULT 1-E PROPERTIES
          ZMATRIX ANGSTROM
          C
          O 1 1.203
          H 1 1.099 2 121.8
          H 1 1.099 2 121.8 3 180.0
          END
          RUNTYPE SCF
          PROPERTY ATOMS
          ENTER

In this example the set of MOs to be used in the property evaluation will be retrieved from that section written in the SCF process, namely section 1 of the Dumpfile i.e. the default section number for the underlying closed-shell SCFTYPE (see Table 1).

A somewhat different approach may be required when computing the one-electron properties derived from a wavefunction with more than one set of MOs (e.g., a UHF wavefunction), or in cases where only the total density matrix, and not an associated set of MOs, is available (e.g., in an MP2 calculation). In both cases, the user may need to ensure that the associated set of spinfree natural orbitals and, where relevant SPIN natural orbitals, are generated by specification of the NATORB directive(s), used to route the NOs to a nominated section on the Dumpfile.

We illustrate this effect by first considering the data requirements when performing a UHF wavefunction. The following data sequence would be required when evaluating the properties based on a direct-UHF calculation, with the computation based on the alpha- and beta-UHF MOs routed to the default sections 1 and 2 respectively under implicit control of the ENTER directive.

          TITLE
          H2CO - 3A2 UHF PROPERTIES - 3-21G BASIS
          MULT 3
          ZMATRIX ANGSTROM
          C
          O 1 1.203
          H 1 1.099 2 121.8
          H 1 1.099 2 121.8 3 180.0
          END
          SCFTYPE DIRECT UHF
          PROPERTY ATOMS
          ENTER

The same calculation may be performed based on the spinfree and spin natural orbitals of the UHF wavefunction; in this case the NATORB data lines will be used to route the spinfree and spin natural orbitals to sections 10 and 11 of the Dumpfile respectively, and these orbitals will be used in computing the 1-electron properties, thus:

          TITLE
          H2CO - 3A2 UHF NO-BASED PROPERTIES - 3-21G BASIS
          MULT 3
          ZMATRIX ANGSTROM
          C
          O 1 1.203
          H 1 1.099 2 121.8
          H 1 1.099 2 121.8 3 180.0
          END
          SCFTYPE DIRECT UHF
          PROPERTY ATOMS
          NATORB 10 
          NATORB SPIN 11
          ENTER

The following data sequence would be required if the user wished to compute the properties of the annihilated UHF wavefunction:

          TITLE
          H2CO - 3A2 annihilated UHF properties 3-21G BASIS 
          MULT 3
          ZMATRIX ANGSTROM
          C
          O 1 1.203
          H 1 1.099 2 121.8
          H 1 1.099 2 121.8 3 180.0
          END
          SCFTYPE DIRECT UHF
          PROPERTY ATOMS
          NATORB 10 ANNIHILATE
          NATORB SPIN 11 ANNIHILATE
          ENTER

Note again that the NOs of the UHF and AUHF wave function are in fact identical, the only difference lying in the occupation numbers.

Now let us consider the date requirements when computing properties at the optimum geometry derived from an MP2 calculation.

          TITLE
          H2CO - X1A1 - MP2 DZ BASIS - PROPERTIES
          ZMATRIX ANGSTROM
          C
          O 1 CO
          H 1 CH 2 HCO
          H 1 CH 2 HCO 3 180.0
          VARIABLES
          CO 1.203\CH 1.099\HCO 121.8
          END
          BASIS DZ
          RUNTYPE OPTIMISE
          PROPERTY ATOMS
          SCFTYPE MP2
          NATORB 20 
          ENTER

Having generated the MP2 optimised geometry, the spinfree natural orbitals will be routed to section 20 on the Dumpfile, and used in the subsequent properties calculation.

13.3 Localised Orbitals

The following data sequence would be required in localising the valence SCF-MOs using the Foster-Boys algorithm, where the LOCAL directive specifies those orbitals deemed to be active in the localisation process.

          RESTART
          TITLE
          H2CO - 3-21G DEFAULT BASIS - VALENCE LMOs
          ZMATRIX ANGSTROM
          C
          O 1 1.203
          H 1 1.099 2 121.8
          H 1 1.099 2 121.8 3 180.0
          END
          RUNTYPE ANALYSE
          LOCAL
          3 TO 8 END
          VECTORS 1
          ENTER 20

Note that the localised orbital module is the only analysis module that creates a new set of eigenvectors, and the user must specify the destination section on the Dumpfile for these orbitals i.e. no default section will be employed. In this case the final set of LMOs will be output to Section 20 of the Dumpfile.

13.4 Distributed Multipole Analysis

The following data sequence would be required in requesting a distributed multipole analysis of the SCF-MOs [32], where the DMA directive instigates the process.

         RESTART
         TITLE
         H2CO - 3-21G DEFAULT BASIS - DMA ANALYSIS
         ZMATRIX ANGSTROM
         C
         O 1 1.203
         H 1 1.099 2 121.8
         H 1 1.099 2 121.8 3 180.0
         END
         RUNTYPE ANALYSE
         DMA
         VECTORS 1
         ENTER

13.5 Graphical Analysis

The following data sequence would be required in generating grids of total density, atom-difference density, electrostatic potential and orbital amplitude for subsequent graphical analysis. The GRAPHICS directive introduces data defining the required graphics processing, with GDEF data defining the grid of points involved, and subsequent CALC and PLOT directives introducing data specifying the required computation associated with the grid (CALC) and corresponding graphical output to be generated (PLOT).

          RESTART
          TITLE
          H2CO - 3-21G DEFAULT BASIS - GRAPHICAL ANALYSIS 
          ZMATRIX ANGSTROM
          C
          O 1 1.203
          H 1 1.099 2 121.8
          H 1 1.099 2 121.8 3 180.0
          END
          RUNTYPE ANALYSE
          GRAPHICS
          GDEF
          TYPE 2D
          POINTS 99
          TITLE
          SQUARE 2D GRID (99*99) 
          CALC
          TYPE ATOM
          TITLE
          H2CO -ATOM DIFFERENCE
          SECTION 150
          PLOT
          TYPE LINE
          TITLE
          ATOM DIFFERENCE DENSITY LINEPRINTER PLOT
          CALC
          TYPE DENS
          SECTION 151
          TITLE
          H2CO - TOTAL DENSITY
          PLOT
          TYPE LINE
          TITLE
          DENSITY LINEPRINTER PLOT
          CALC
          TYPE MO 2
          TITLE
          H2CO MO 2 AMPLITUDE
          SECTION 152
          PLOT
          TYPE LINE
          TITLE
          MO 2 LINEPRINTER PLOT
          GDEF
          TYPE 2D
          POINTS 25
          TITLE 
          SQUARE 2D GRID (25*25)
          CALC
          TYPE POTE
          TITLE
          H2CO - POTENTIAL
          SECTION 153
          PLOT
          TYPE LINE
          TITLE
          POTENTIAL LINEPRINTER PLOT 
          VECTORS 1
          ENTER

The resolution of each plot is controlled by the size of the grid, via the POINTS sub-directive of GDEF. Note that the TYPE sub-directive of CALC defines the type of grid (ATOM, DENS, MO and POTE for atom-difference, total density, orbital amplitude and electrostatic potential respectively). In the present example output is restricted to the line-printer, through the LINE parameter in the PLOT data.

13.6 Population Analysis

The following data sequence would be required in performing an extended population analysis of the valence SCF-MOs, where the MULLIKEN directive specifies those orbitals for which printed output is required. The ATOM and ORBITAL keyword request the emphasis in the analysis generated through the grouping of basis functions [34]. The sequence of integers specified on the MULLIKEN line specifies those MOs for which printed output is required.

          RESTART
          TITLE
          H2CO - 3-21G BASIS - ANALYSIS OF VALENCE MOs
          ZMATRIX ANGSTROM
          C
          O 1 1.203
          H 1 1.099 2 121.8
          H 1 1.099 2 121.8 3 180.0
          END
          RUNTYPE ANALYSE
          MULLIKEN ATOM ORBITAL 3 TO 8 END
          VECTORS 1
          ENTER

Note that it is also possible to define the groups of basis functions through user input. Thus the following data would perform the same analysis as the ATOM specification above, where the GROUP keyword on the MULLIK data line indicates that subsequent data lines will follow, terminated by the END keyword, that will assign the basis functions to user-defined groups.

          RESTART
          TITLE
          H2CO - 3-21G BASIS - INPUT GROUPS FOR ANALYSIS
          ZMATRIX ANGSTROM
          C
          O 1 1.203
          H 1 1.099 2 121.8
          H 1 1.099 2 121.8 3 180.0
          END
          RUNTYPE ANALYSE
          MULLIKEN GROUP  3 TO 8 END
          CATOM 1 TO 9
          OATOM 10 TO 18
          H1ATOM 19 20
          H2ATOM 21 22
          END
          VECTORS 1
          ENTER

13.7 Morokuma Energy Decomposition Analysis

The following example illustrates how the Morokuma EDA is performed using a sequence of three separate GAMESS-UK input decks, one for each of the two fragments and an analysis job for the supermolecule.

The Class 2 MOROKUMA directive controls the job, and may take one of two forms, depending on whether a fragment SCF or an interaction calculations required. The directive sequence:

          MOROKUMA FRAG NUMBER TAG

specifies the RHF calculation on one of the fragments. NUMBER should be 1 or 2, indicating the position of the fragment in the supermolecule. TAG is replaced with a string to identify the fragment; the job will result in a file of this name containing the fragment basis and wavefunction information being written in the working directory of the job. There is currently an 8 character limit on TAG. The sequence:

          MOROKUMA INTERACT TAG1 TAG2

requests that an interaction energy analysis be performed. The geometry is assumed to be that of the supermolecule, and the two tags denote the fragment files from two previous runs under control of MOROKUMA FRAG as above.

A number of restrictions should be noted when using the morokuma analysis module:

• The implementation is restricted to RHF calculations using the conventional (non-direct) SCF module.
• The use of symmetry (including symmetry adaption) must be disabled for all component jobs. The atoms in the supermolecule must be presented in the same order as that obtained by concatenating the two fragments and the basis sets specified for the separate tasks must correspond.
• The code is developmental, although it is believed to work correctly within the above limits, prospective users are advised to contact the authors.
• Morokuma EDA jobs cannot be restarted.

          TITLE 
          MOROKUMA TEST FRAG2
          ADAPT OFF
          NOSYM
          GEOMETRY
           0.00000000  -1.10092542 -1.43475395  1.0 H
           0.00000000  -1.10092542  1.43475395  1.0 H
           0.00000000   0.00000000  0.00000000  8.0 O
          END
          MOROKUMA FRAG 1 FRAG1
          BASIS SV 4-31G
          ENTER


          TITLE 
          MOROKUMA TEST
          ADAPT OFF
          NOSYM
          GEOMETRY
           3.24201636   2.02583666  0.00000000  1.0 H
           4.24693920   4.71362490  0.00000000  1.0 H
           4.77568401   2.98417857  0.00000000  8.0 O
          END
          MOROKUMA FRAG 2 FRAG2
          BASIS SV 4-31G
          ENTER


          TITLE 
          MOROKUMA TEST
          ADAPT OFF
          NOSYM
          GEOMETRY
           0.00000000  -1.10092542 -1.43475395  1.0 H
           0.00000000  -1.10092542  1.43475395  1.0 H
           0.00000000   0.00000000  0.00000000  8.0 O
           3.24201636   2.02583666  0.00000000  1.0 H
           4.24693920   4.71362490  0.00000000  1.0 H
           4.77568401   2.98417857  0.00000000  8.0 O
          END
          BASIS SV 4-31G
          MOROK INTERACT FRAG1 FRAG2
          VECTORS ATOMS
          ENTER