IOp(3/110=…)

Supplementary User’s Guide

for a version of GAUSSIAN03

(rev. D.01)

modified to allow nuclear charge stabilization calculations for metastable states

This version of GAUSSIAN03

has been modified and tested by

Prof. Piotr Skurski[1] while working in

The Simons Group at The

Henry Eyring Center for Theoretical Chemistry

Department of Chemistry, University of Utah

Salt Lake City, UT 84112, U.S.A.

Introduction

When studying metastable electronic states that can decay by ejection of an electron, one cannot straightforwardly employ variation-based electronic structure tools[1]. For such cases, the lowest-energy state corresponds to that of a free electron (infinitely distant and with zero velocity) plus a system with one fewer electrons. Standard methods thus suffer so-called variational collapse and tend to converge to such “free-electron” descriptions. Let’s consider an example offered by the formamide molecule near its equilibrium geometry. The p-symmetry molecular orbitals (MOs) are depicted qualitatively below where their nodal characteristics are emphasized. The lowest two p MOs describe the delocalized p bonding and non-bonding orbitals. The unoccupied MO is the anti-bonding p* orbital.

An SCF calculation using a correlation consistent augmented polarized valence double-zeta (aug-cc-pVDZ) basis set produces the orbitals shown below. The orbital energies for the bonding and non-bonding p MOs (labeled HOMO-2 and HOMO) are -15.4 and -11.5 eV, respectively. The HOMO-1 orbital is a lone pair orbital on the oxygen atom.

The SCF orbital energy of the lowest unoccupied molecular orbital (LUMO) is 0.72 eV, which suggests that an electron of 0.72 eV might attach to produce the formamide anion. However, as the picture shown below illustrates, the LUMO is not even of p* symmetry, nor is the LUMO + 1 or the LUMO + 2 orbital. In fact, these three unoccupied orbitals do not have any significant valence character; most of their amplitude is outside the formamide molecule’s molecular skeleton. They are, within the finite atomic orbital basis used, approximations to the free-electron orbital.

The lowest unoccupied orbital of p* character is the LUMO + 3, which is also shown below, and this orbital has an energy of + 2.6 eV. However, in a different atomic orbital basis, the lowest unoccupied orbital of p* symmetry would not necessarily be LUMO + 3, nor would it necessarily have an orbital energy near 2.6 eV.

To give the reader an idea of how these orbitals look for other atomic orbital basis sets, we show below four of the orbitals obtained when a 6-31G** basis is employed.

We see that the HOMO-2 and HOMO are still p bonding and non-bonding orbitals and HOMO-1 is still an oxygen lone pair orbital, but now LUMO is the p* anti-bonding orbital. It is important to notice that the desired p* orbital may be the LUMO in one basis but might be another orbital in a different basis as it is in the examples shown above.

Clearly, estimating the vertical electron affinity (EA) by using the energy of the LUMO is wrong, but we need to make clear what is wrong and how to properly identify the correct orbital. In fact, most of the low-energy vacant orbitals are nothing but attempts, within the finite orbital basis used, to represent a free electron plus a neutral formamide molecule. This illustrates the variational collapse problem mentioned above.

Let us now explain how we can fix this error. To focus on the problem, we consider the effective potential that an electron approaching a formamide molecule from afar would experience. Also keep in mind that, as is conventional in electron scattering calculations, we may wish to decompose the wave function describing this “attached” electron into products of radial and angular terms:

Y(r,q,f) = S_L,My_L,M(r) Y_L,M(q,f).

Substituting such an expansion into the Schrödinger equation

{- h²/2m r ^-2∂/∂r(r²∂/∂r) + h²L²/2mr² + V} Y = E Y

produces the usual set of coupled-channel equations. Multiplying the Schrödinger equation on the left by Y*_L,M and integrating over the angles q and f, one obtains a radial equation for the y_L,M component of the wave function:

{- h²/2m r^-2∂/∂r(r²∂/∂r) + h²<L²>/2mr² + <V>}y_L,M = Ey_L,M

where <V> denotes the angular average of the electron-molecule potential

<V> = ò Y*_L,MV Y_L,Msinq dq df.

When the distance r of the electron from the molecule is large, <V> is dominated by an attractive Coulomb factor if the electron is interacting with a cation, by a charge-dipole potential if the electron is interacting with a polar neutral molecule, and by a repulsive Coulomb factor if the electron is interacting with a negative ion. However, the centrifugal potential h²<L²>/2mr²always varies as r^-2.

For the formamide case at hand, the longest-range component of <V> is the charge-dipole term, which varies as r^-2. Moreover, the nodal character of the p* orbital into which the incoming electron is to attach has dominant L = 3 character. To see this, we view this orbital from a long distance at which its three nuclear centers are nearly on top of one another as shown below.

When viewed as having the O, C, and N nuclei on top of one another, this p* orbital clearly has nodal properties like that of an f-orbital which is why L = 3 is dominant at large-r.

An electron in an orbital having angular momentum L experiences an effective radial potential (i.e., the sum of <V> and the centrifugal potential) that varies as shown qualitatively below for a neutral molecule. For an electron interacting with an anion, the repulsive long-range part of the potential would also include a Coulomb term e²/r. In such cases, the barrier that acts to constrain the electron in the metastable state arises from both Coulomb and centrifugal factors. Nevertheless, the methods discussed here can still be used to characterize the metastable state’s energy and wave function.

If, as in the case of formamide, the component of the potential generated by the attractive valence-range influences of the O, C, and N centers is not strong enough, no bound state of the anion will exist[2]. For such a case, only the metastable so-called L = 3 shape resonance state will occur and it will have an energy (the heavy horizontal line) and a radial wave function as shown below.

In the valence region, the resonance function has large amplitude, suggesting that the electron is rather localized, it decays exponentially in the classically forbidden tunneling region, and it has sinusoidal oscillations beyond this region with the local de Broglie wavelength relating to the electron’s kinetic energy.

However, at energies both above and below that of the shape resonance state, there exists a continuum of other states. Those lying below the resonance energy vary with r as shown below.

They have little amplitude in the valence region and have large amplitude at large-r and they have longer de Broglie wavelength and thus lower energy than the resonance state. There are also non-resonant states lying energetically above the shape resonance; they have little amplitude in the valence region and larger amplitude at large-r but with shorter de Broglie wavelength corresponding to higher kinetic energy.

So, how can we identify and characterize the resonance wave function (and its energy) when it is “buried” within a continuum of higher- and lower- energy states? The example discussed earlier makes it clear we cannot rely on the orbital energy nor on choosing the LUMO to find the desired state. In the nuclear charge stabilization method[3], we attempt to smoothly slightly enhance the valence-region attractive character of the potential V that the attached electron experiences to an extent that pulls the energy of the metastable resonance state below zero thus rendering it stable[4]. We do this by altering the nuclear charges of those atoms over which the resonance MO is expected to be localized in the valence region. For formamide, we alter the O, C, and N nuclear charges from their actual values of 8, 6, and 7, respectively, to, for example, 8+dq, 6+dq, and 7+dq, respectively. The nuclear charges appear in the electronic Hamiltonian in the electron-nuclear interaction potential

V_e,n= S_j=1,NS_K=1.M(-Z_Ke²/|r_j-R_K|)

where the sum over j runs over all N electrons, the sum over K runs over all M nuclei of charge Z_K. So, the nuclear charge scaling in formamide introduces an additional one-electron potential

dV_e,n= - dqS_j=1,N[(e²/|r_j-R_O|) (e²/|r_j-R_C|) + (e²/|r_j-R_N|)

that acts to differentially stabilize the electron’s potential energy near these nuclear centers. Visually, we can represent the effects of the scaled nuclear charges on the electron-molecule potential as shown below.

Here, in black, we reproduce the original unscaled potential and its resonance wave function. In blue and red, we illustrate potentials in which the nuclear charges have been incremented by small (blue) or somewhat larger (red) amounts. The key thing to notice that is, if the nuclear charge enhancement is large enough, the valence-range component of the potential will be lowered enough to render the resonance state bound (notice the lack of continuum components at large-r and notice that the energy lies below zero) rather than metastable and thus amenable to studying using conventional variational-based quantum chemistry tools. In this manner, one can carry out conventional calculations on the nuclear charge enhanced species for a series of dq values (all of which must be large enough to render the desired state bound) and then extrapolate to dq ® 0 to allow the resonance state to be identified from the finite-dq calculations’ data.

In practice, to use the nuclear charge stabilization method, one (for a series of dq values)

1. Identifies those nuclei over which the valence component of the desired resonance state’s orbital will be distributed;

2. Modifies the nuclear charges of these nuclei (one need not use equal charge increments for all the nuclei[5]) by increasing them by an amount dq;

3. Carries out a standard, bound-state, quantum calculation (SCF, MPn, coupled-cluster, or whatever) on the electron-attached[6] and non-attached states using the scaled nuclear charges to obtain attached (E*) and non-attached (E) state energies, after which

4. One plots the energy difference (E-E*) vs. dq but using only dq values large enough to render E-E* > 0 (i.e., to make the attached state bound relative to the non-attached state), and

5. One then extrapolates the plot to dq ® 0 to obtain an estimate of the energy of the electron-attached state relative to that of the parent non-attached species (i.e., in the extrapolation, one will find E-E* negative, meaning the electron attached species lies above its parent).

In this manner, one achieves a nuclear charge stabilization estimate of the energy of the desired metastable state.

In the figure shown below, we show Koopmans’ theorem[7], SCF-level, MP2-level, and coupled-cluster level energy differences (D) for the formamide and formamide p* anion obtained using the aug-cc-pVDZ basis set discussed earlier. Only for values of dq equal to 0.18 or larger was the anion found to be electronically stable, so only data obtained with dq = 0.18 and higher were used in constructing these plots.

The dq ® 0 extrapolated energies for each case are as follows:

D = –4.3 eV (KT, red) and the fit’s correlation coefficient (squared) is: r²=0.99957.

D = –3.6 eV (SCF, blue), r²=0.99996

D = –3.1 eV (MP2, green), r²=0.99994

D = –3.1 eV (CCSD(T), magenta), r²=0.99994.

When analogous nuclear charge stabilization calculations are performed using the 6-31G** basis, the plot shown below results.

For this basis set, the extrapolated predictions of the D-values are as follows:

D(KT) = -5.6 eV

D(SCF) = -4.3 eV

D(MP2) = -4.4 eV

D(CCSD(T)) = -4.4 eV.

These extrapolations represent the nuclear charge stabilization’s predictions for the energy of the anion’s metastable p* resonance state relative to that of the neutral for each level of theory.

In addition, below we show the p* LUMO obtained for the aug-cc-pVDZ basis set for dq = 0.19 and 0.26 to show how it differs significantly from the LUMO obtained in the non-scaled calculation which, as we explained earlier, cannot be trusted to relate to the desired resonance state.

Clearly, this orbital has the desired valence character and nodal pattern. For dq = 0.26 it is more compact than for dq = 0.19 because of the enhanced nuclear attraction in the former case. For the 6-31G** basis, this same orbital is shown below also for dq = 0.16 and dq = 0.29.

We see that the qualitative character of the desired p* orbital does not depend on the basis set employed although the quantitative values of the resonance state energy do. For the larger basis and using the highest-level of theory (the CCSD(T) data), the p* attached anion is predicted to lie 3.1 eV above the neutral formamide at the neutral’s equilibrium geometry. Experiments using so-called electron transmission spectroscopy methods find[8] a resonance state to lie ca. 2 eV above the neutral, so one can see that obtaining accurate estimates of the energies of such metastable states is difficult even when rather good atomic basis sets are employed.

Details for using the G03 modifications.

There are two steps that must be taken; one must identify those nuclei whose charges are to be altered and one must tell by how much to increase or decrease the charge of each nucleus.

Specifiying the nuclear centers whose charges are to be modified

IOp(3/110=…) and IOp(7/110=…)

These IOps are used to tell Gaussian the highest numbered nuclear center range of atoms (centers in the geometry input list in your Gaussian .com file) in which there are the nuclei whose charges are to be modified.

Example #1:

To modify only the charge on center number 1, one needs to specify IOp(3/110=1) and IOp(7/110=1).

Example #2:

To modify two nuclei whose numbers (in the geometry input specification) are 1 and 2, one needs to specify IOp(3/110=2) and IOp(7/110=2).

Example #3:

To modify four nuclei whose numbers (in the geometry input specification) are 13, 26, 27, and 67, one needs to specify IOp(3/110=67) and IOp(7/110=67).

Briefly, the values of option 110 in overlay 3 and option 110 in overlay 7 have to be set to the highest numbered center that is to be modified, no matter how many centers will actually be modified. Note that option 110 should have the same value in both overlays (i.e., IOp(3/110=X₁) and IOp(7/110=X₂), where X₁=X₂).

Warning: the largest possible value for this option is 90 in both overlays (see the “Limitations and known problems” section at the end of this manual), so make sure you number your atoms so that all of the atoms you intend to subject to charge scaling occur within the first 90 atoms.

Modifying the nuclear charges on the selected centers

IOp(3/111…200) and IOp(7/111…200)

These IOps are used both to denote the centers to be modified and to specify the additional charge (dq) that is to be added to that center (or subtracted from it).

To modify nuclear charge on center #N, the IOptions 3/(110+N) and 7/(110+N) must be set to a proper value. The additional charge to be added to (or subtracted from) the N^th nucleus must be multiplied by 10⁵ (100000) and set as the value in both IOptions 3/(110+N) and 7/(110+N).

Example #1:

To add the charge of 0.1 a.u. to the center number 1, one needs to specify IOp(3/111=10000) and IOp(7/111=10000) because:

ü the center number is 1 which corresponds to options 111 (110+1) in overlays 3 and 7 , and

ü the charge to be added is 0.1 a.u. which corresponds to values of 10000 for options 111 in overlays 3 and 7 (because 0.1×10⁵ = 10000).

Example #2:

To modify two nuclei whose numbers in the geometry input specification are 4 and 17, by adding an extra positive charge of 0.03 a.u. to nucleus on center 4 and charge of 0.58 a.u. to the nucleus on center 17, one needs to specify IOp(3/114=3000,3/127=58000) IOp(7/114=3000,7/127=58000) because

ü the number of the first center to be modified is 4 which corresponds to options 114 (110+4) in overlays 3 and 7 , and

ü the number of the second center to be modified is 17 which corresponds to options 127 (110+17) in overlays 3 and 7 , and

ü the charge to be added to the first modified center is 0.03 a.u. which corresponds to values of 3000 for options 114 in overlays 3 and 7 (because 0.03×10⁵ = 3000) for this center, and

ü the charge to be added to the second modified center is 0.58 a.u. which corresponds to values of 58000 for options 127 in overlays 3 and 7 (because 0.58×10⁵ = 58000) for this center.

Recall, that even though you are choosing the center numbers here, you still need to set the proper range by using IOp(3/110=…) and IOp(7/110=…)

(see the preceding section).

Input examples

Here we present a few examples of Gaussian03 input files showing single-point energy calculations and geometry optimization jobs that involve modified nuclear charges.

Example #1:

%chk=testSP

#p RHF 6-31G(d) SCF=(Tight,NoVarAcc,Maxcycle=512,IntRep)

GFInput IOp(6/7=3) IOp(3/32=1)

IOp(3/110=3,3/111=33000,3/112=34000,3/113=33000)

IOp(7/110=3,7/111=33000,7/112=34000,7/113=33000)

Test job #1

0 1

C 1 co2

N 2 nc3 1 nco3

H 2 hc4 1 hco4 3 dih4

H 3 hn5 2 hnc5 1 dih5

H 3 hn6 2 hnc6 1 dih6

co2 1.277811

nc3 1.496812

nco3 116.154

hc4 1.002000

hco4 135.119

dih4 172.107

hn5 1.004634

hnc5 107.935

dih5 -119.695

hn6 1.002000

hnc6 112.357

dih6 0.000

In this job, the first three atoms are modified: charges of 0.33, 0.34, and 0.33 a.u. are added to the oxygen atom, carbon atom, and nitrogen atom, respectively. So the resulting charges for atoms O, C, N are 8.33, 6.34, and 7.33 a.u., respectively, while those for hydrogen atoms remain unchanged. This is achieved in the 4^th and 5^th line of the input file. The Hartree-Fock single point energy for such a system will be calculated, with the 6-31G(d) basis sets, for the given geometry.

Example #2:

%chk=testOPT

#p UHF 6-31G(d) SCF=(Tight,NoVarAcc,Maxcycle=512,IntRep)

GFInput IOp(6/7=3) IOp(3/32=1) OPT GUESS=READ

IOp(3/110=3,3/111=33000,3/112=34000,3/113=33000)

IOp(7/110=3,7/111=33000,7/112=34000,7/113=33000)

Test job #2

-1 2

C 1 co2

N 2 nc3 1 nco3

H 2 hc4 1 hco4 3 dih4

H 3 hn5 2 hnc5 1 dih5

H 3 hn6 2 hnc6 1 dih6

co2 1.277811

nc3 1.496812

nco3 116.154

hc4 1.002000

hco4 135.119

dih4 172.107

hn5 1.004634

hnc5 107.935

dih5 -119.695

hn6 1.002000

hnc6 112.357

dih6 0.000

In this job, the first three atoms are modified: charges of 0.33, 0.34, and 0.33 a.u. are added to the oxygen atom, carbon atom, and nitrogen atom, respectively. So the resulting charges for atoms O, C, N are 8.33, 6.34, and 7.33 a.u., respectively, while those for hydrogen atoms remain unchanged. This is achieved in the 4^th and 5^th line of the input file. The Hartree-Fock geometry optimization for such an anion will be performed, with the 6-31G(d) basis sets.

Example #3:

%chk=testOPT2

#p UHF aug-cc-pVDZ SCF=(Tight,NoVarAcc,Maxcycle=512)

GFInput IOp(6/7=3) IOp(3/32=1) OPT GUESS=READ

IOp(3/110=5,3/111=33000,3/112=34000,3/115=33000)

IOp(7/110=5,7/111=33000,7/112=34000,7/115=33000)

Test job #3

-1 2

O 0.000000 0.000000 0.000000

C 0.000000 0.000000 1.277811

H -0.700350 0.097094 1.987802

H 1.412538 -0.830293 2.498962

N 1.343557 0.000000 1.937584

H 2.094138 0.000000 1.273783

In this job, there are also three modified atoms: charges of 0.33, 0.34, and 0.33 a.u. are added to the oxygen atom, carbon atom, and nitrogen atom, respectively. So the resulting charges for atoms O, C, N are 8.33, 6.34, and 7.33 a.u., respectively, while those for hydrogen atoms remain unchanged. This is achieved again in the 4^th and 5^th line of the input file. The Hartree-Fock geometry optimization for such an anion will be performed, with the 6-31G(d) basis sets. (Please note how the centers for modifications are specified via IOps and note that the atoms whose charges have been changed are the first, second, and fifth atoms).

Some limitations and known problems

The maximum number of centers whose charges can be modified in one job is limited to 90, which is limited by the array size used for storing the IOps.

The maximum number of atoms in a molecule is not limited other than being subject to the same limitations as in any regular Gaussian03 input files.

The highest possible number for a modified center is 90. This means that the centers must be numbered so that all those to be charge-modified occur in the first 90.

The modified version of Gaussian03 has been tested for single-point energy and geometry optimization jobs. Calculating harmonic vibrational frequencies is also possible. Note that changing the nuclear charges does not affect the atomic masses, so quantities that depend on atomic masses (such as vibrational frequencies) are calculated in a standard G03 manner.

The value of the nuclear-nuclear repulsion energy, as well as its contributions to the gradients and Hessians, are calculated using the conventional (non-modified) nuclear charges.

Some modifications of the nuclear charges are to be avoided because they can produce unexpected errors. In particular, changing any nuclear charge by 1 a.u. (or more) leads to problems since, de facto, this changes the chemical identity of the element. In such cases errors reported by MinBas appear (due to the comparison made by Gaussian03 in which the atomic number is compared with its nuclear charge).

In some cases, both the SCF procedure and the geometry optimization process can be more difficult to converge when modified nuclear charges are used, so it is recommended that one use the MaxCycle keyword to increase the number of maximum SCF iterations and optimization steps to improve chances for successful convergence.

Molecular properties such as populations (e.g., Mulliken population analysis, ESP-fitted charges according to the Merz-Singh-Kollman scheme, etc.), multipole moments, polarizability tensors, and other quantities that depend on the nuclear charges are evaluated with the modified nuclear charges.

Please report any problems to:

Piotr Skurski (piotr@hydrogen.hec.utah.edu).

[1] Chemistry Department, University of Gdansk, Gdansk, Poland