# Category Archives: White papers

## Population Analysis in the Excited State with Gaussian

To calculate what the bonding properties of a molecule are in a particular excited state we can run any population analysis following the root of interest. This straightforward procedure takes two consecutive calculations since you don’t necessarily know before hand which excited state is the one of interest.

The regular Time Dependent Density Functional Theory (TD-DFT) calculation input with Gaussian 16 looks as follows (G09 works pretty much the same), let us assume we’ve already optimized the geometry of a given molecule:

```%OldChk=filename.chk
%nprocshared=16
%chk=filename_ES.chk

#p TD(NStates=10,singlets) wb97xd/cc-pvtz geom=check guess=read

Title Card Required

0 1
--blank line--```

This input file retrieves the geometry and wavefunction from a previous calculation from filename.chk and doesn’t write anything new into it (that is what %OldChk=filename.chk means) and creates a new checkpoint where the excited states are calculated (%chk=filename_ES.chk)

In the output you search for the transition which peeks your interest; most often than not you’ll be interested in the one with the highest oscillator strength, f. The oscillator strength is a dimensionless number that represents the ratio of the observed, integrated, absorption coefficient to that calculated for a single electron in a three-dimensional harmonic potential [Harris & Bertolucci, Symmetry and Spectroscopy]; in other words, it is related to the probability of that transition to occur, and therefore it takes values from 0.0 to 1.0 (for single photon absorption processes.)

The output of this calculation looks as follows, the value of f for every excitation is reported together with its energy and the orbital transitions which comprise it.

``` Excitation energies and oscillator strengths:

Excited State   1:      Singlet-A      3.1085 eV  398.86 nm  f=0.0043  <S**2>=0.000
56 -> 59        -0.11230
58 -> 59         0.69339
This state for optimization and/or second-order correction.
Total Energy, E(TD-HF/TD-DFT) =  -1187.56377917
Copying the excited state density for this state as the 1-particle RhoCI density.

Excited State   2:      Singlet-A      4.0827 eV  303.68 nm  f=0.0016  <S**2>=0.000
52 -> 59         0.46689
52 -> 64        -0.20488
53 -> 59         0.19693
54 -> 59         0.40414
54 -> 64        -0.16261
...
...
Excited State   8:      Singlet-A      5.2345 eV  236.86 nm  f=0.8063  <S**2>=0.000
52 -> 60         0.17162
53 -> 59         0.47226
53 -> 60        -0.11771
54 -> 59        -0.27658
54 -> 60        -0.22006
55 -> 59         0.20496
56 -> 59         0.15029
```

Now we’ve selected excited state #8 because it has the largest value of f from the lot, we use the following input to read in the geometry from the old checkpoint file and we generate a new one in case we need it for something else. The input file for doing all this looks as follows (I’ve selected as usual the Natural Bond Orbital population analysis):

```%oldchk=a_ES.chk
%nprocshared=16
%chk=a_nbo.chk

Title Card Required

0 1

\$NBO BOAO BNDIDX E2PERT \$END

--blank line--

```

The flags at the bottom request the calculation of Wiberg Bond Indexes (BNDIDX) as well as Bond Order in the Atomic Orbital basis (BOAO) and a second order perturbation theory for the electronic delocalization (E2PERT). Now we can compare the population analysis between ground and the 8th excited state; check figure 1 and notice the differences in Wiberg’s bond order for this complex made of two molecules and one Na+ cation. Figure 1. Natural Population Analysis comparison for a supramolecular arrangement. Numbers next in brackets correspond to the sum of charges for each molecule. Notice the significant change in charges for each molecule when going from the ground to the 8th excited state.

In this example we can observe that in the ground state we have a neutral and a negative molecule together with a Na+ cation, but when we analyze the population in the 8th excited state both molecules acquire a similar charge, ca. 0.46e, which means that some of the electron density has been transferred from the negative one to the neutral molecule, forming an Electron Donor-Acceptor complex (EDA) in the excited state.

This procedure can be extended to any other kind of population analysis and their derived combination, e.g. one could calculate their condensed fukui functions in the Nth excited state; but beware! These calculations yield vertical excitations, should the excited state of interest have a minimum we can first optimize the ES geometry and then perform the population analysis on said geometry; just add the opt keyword to perform both jobs in one go, but bear in mind that the NBO population analysis is performed before and after the optimization process so look for the tables and values closer to the end of the output file.

In the case of open shell systems the procedure is the same but one should be extremely careful in searching for the total population analysis since the output file contains this table for the alpha and beta populations separately as well as the added values for the total number of electrons.

## Percentage of Molecular Orbital Composition – G09,G16

Canonical Molecular Orbitals are–by construction–delocalized over the various atoms making up a molecule. In some contexts it is important to know how much of any given orbital is made up by a particular atom or group of atoms, and while you could calculate it by hand given the coefficients of each MO in terms of every AO (or basis set function) centered on each atom there is a straightforward way to do it in Gaussian.

If we’re talking about ‘dividing’ a molecular orbital into atomic components, we’re most definitely talking about population analysis calculations, so we’ll resort to the pop keyword and the orbitals option in the standard syntax:

`#p M052x/cc-pVDZ pop=orbitals`

This will produce the following output right after the Mulliken population analysis section:

```Atomic contributions to Alpha molecular orbitals:
Alpha occ 140 OE=-0.314 is Pt1-d=0.23 C38-p=0.16 C31-p=0.16 C36-p=0.16 C33-p=0.15
Alpha occ 141 OE=-0.313 is Pt1-d=0.41
Alpha occ 142 OE=-0.308 is Cl2-p=0.25
Alpha occ 143 OE=-0.302 is Cl2-p=0.72 Pt1-d=0.18
Alpha occ 144 OE=-0.299 is Cl2-p=0.11
Alpha occ 145 OE=-0.298 is C65-p=0.11 C58-p=0.11 C35-p=0.11 C30-p=0.11
Alpha occ 146 OE=-0.293 is C58-p=0.10
Alpha occ 147 OE=-0.291 is C22-p=0.09
Alpha occ 148 OE=-0.273 is Pt1-d=0.18 C11-p=0.12 C7-p=0.11
Alpha occ 149 OE=-0.273 is Pt1-d=0.18
Alpha vir 150 OE=-0.042 is C9-p=0.18 C13-p=0.18
Alpha vir 151 OE=-0.028 is C7-p=0.25 C16-p=0.11 C44-p=0.11
Alpha vir 152 OE=0.017 is Pt1-p=0.10
Alpha vir 153 OE=0.021 is C36-p=0.15 C31-p=0.14 C63-p=0.12 C59-p=0.12 C38-p=0.11 C33-p=0.11
Alpha vir 154 OE=0.023 is C36-p=0.13 C31-p=0.13 C63-p=0.11 C59-p=0.11
Alpha vir 155 OE=0.027 is C65-p=0.11 C58-p=0.10
Alpha vir 156 OE=0.029 is C35-p=0.14 C30-p=0.14 C65-p=0.12 C58-p=0.11
Alpha vir 157 OE=0.032 is C52-p=0.09
Alpha vir 158 OE=0.040 is C50-p=0.14 C22-p=0.13 C45-p=0.12 C17-p=0.11
Alpha vir 159 OE=0.044 is C20-p=0.15 C48-p=0.14 C26-p=0.12 C54-p=0.11
```

Alpha and Beta densities are listed separately only in unrestricted calculations, otherwise only the first is printed. Each orbital is listed sequentially (occ = occupied; vir = virtual) with their energy value (OE = orbital energy) in atomic units following and then the fraction with which each atom contributes to each MO.

By default only the ten highest occupied orbitals and ten lowest virtual orbitals will be assessed, but the number of MOs to be analyzed can be modified with orbitals=N, if you want to have all orbitals analyzed then use the option AllOrbitals instead of just orbitals. Also, the threshold used for printing the composition is set to 10% but it can be modified with the option ThreshOrbitals=N, for the same compound as before here’s the output lines for HOMO and LUMO (MOs 149, 150) with ThreshOrbitals set to N=1, i.e. 1% as occupation threshold (ThreshOrbitals=1):

```Alpha occ 149 OE=-0.273 is Pt1-d=0.18 N4-p=0.08 N6-p=0.08 C20-p=0.06 C13-p=0.06 C48-p=0.06 C9-p=0.06 C24-p=0.05 C52-p=0.05 C16-p=0.04 C44-p=0.04 C8-p=0.03 C15-p=0.03 C17-p=0.03 C45-p=0.02 C46-p=0.02 C18-p=0.02 C26-p=0.02 C54-p=0.02 N5-p=0.01 N3-p=0.01
Alpha vir 150 OE=-0.042 is C9-p=0.18 C13-p=0.18 C44-p=0.08 C16-p=0.08 C15-p=0.06 C8-p=0.06 N6-p=0.04 N4-p=0.04 C52-p=0.04 C24-p=0.04 N5-p=0.03 N3-p=0.03 C46-p=0.03 C18-p=0.03 C48-p=0.02 C20-p=0.02```

The fragment=n label in the coordinates can be used as in BSSE Counterpoise calculations and the output will show the orbital composition by fragments with the label "Fr", grouping all contributions to the MO by the AOs centered on the atoms in that fragment.

As always, thanks for reading, sharing, and rating. I hope someone finds this useful.

## Fixing the error: Bad data into FinFrg

I found this error in the calculation of two interacting fragments, both with unpaired electrons. So, two radicals interact at a certain distance and the full system is deemed as a singlet, therefore the unpaired electron on each fragment have opposite spins. The problem came when trying to calculate the Basis Set Superposition Error (BSSE) because in the Counterpoise method you need to assign a charge and multiplicity to each fragment, however it’s not obvious how to assign opposite spins.

The core of the problem is related to the guess construction; normally a Counterpoise calculation would look like the following example:

```#p B3LYP/6-31G(d,p) counterpoise=2

-2,1 -1,2 -1,2
C(Fragment=1)        0.00   0.00   0.00
O(Fragment=2)        1.00   1.00   1.00
...```

In which the first pair of charge-multiplicity numbers correspond to the whole molecule and the following to those of each fragment in increasing order of N (in this case, N = 2). So for this hypothetical example we have two anions (but could easily be two cations) each with an unpaired electron, yielding a complex of charge = -2 and a singlet multiplicity which implies those two unpaired electrons have opposite spin. But if the guess (the initial trial wavefunction from which the SCF will begin) has a problem understanding this then the title error shows up:

```Bad data into FinFrg
Error termination via Lnk1e ...```

The solution to this problem is as simple as it may be obscure: Create a convenient guess wavefunction by placing a negative sign to the multiplicity of one of the fragments in the following example. You may then use the guess as the starting point of other calculations since it will be stored in the checkpoint file. By using this negative sign we’re not requesting a negative multiplicity, but a given multiplicity of opposite spin to the other fragment.

```#p B3LYP/6-31G(d,p) guess=(only,fragment=2)

-2,1 -1,2 -1,-2
C(Fragment=1)        0.00   0.00   0.00
O(Fragment=2)        1.00   1.00   1.00
...```

This way, the second fragment will have the opposite spin (but the same multiplicity) as the first fragment. The only keyword tells gaussian to only calculate the guess wave function and then exit the program. You may then use that guess as the starting point for other calculations such as my failed Counterpoise one.

## Geometry Optimizations for Excited States

Electronic excitations are calculated vertically according to the Frank—Condon principle, this means that the geometry does not change upon the excitation and we merely calculate the energy required to reach the next electronic state. But for some instances, say calculating not only the absorption spectra but also the emission, it is important to know what the geometry minimum of this final state looks like, or if it even exists at all (Figure 1). Optimizing the geometry of a given excited state requires the prior calculation of the vertical excitations whether via a multireference method, quantum Monte Carlo, or the Time Dependent Density Functional Theory, TD-DFT, which due to its lower computational cost is the most widespread method.

Most single-reference treatments, ab initio or density based, yield good agreement with experiments for lower states, but not so much for the higher excitations or process that involve the excitation of two electrons. Of course, an appropriate selection of the method ensures the accuracy of the obtained results, and the more states are considered, the better their description although it becomes more computationally demanding in turn. Figure 1. The vertical excitation does not match the minimum on the excited state

In Gaussian 09 and 16, the argument to the ROOT keyword selects a given excited state to be optimized. In the following example, five excited states are calculated and the optimization is requested upon the second excited state. If no ROOT is specified, then the optimization would be carried out by default on the first excited state (Where L.O.T. stands for Level of Theory).

`#p opt TD=(nstates=5,root=2) L.O.T.`

Gaussian16 includes now the calculation of analytic second derivatives which allows for the calculation of vibrational frequencies for IR and Raman spectra, as well as transition state optimization and IRC calculations in excited states opening thus an entire avenue for the computation of photochemistry.

If you already computed the excited states and just want to optimize one of them from a previous calculation, you can read the previous results with the following input :

`#p opt TD=(Read,Root=N) L.O.T. Density=Current Guess=Read Geom=AllCheck`

Common problems. The following error message is commonly observed in excited state calculations whether in TD-DFT, CIS or other methods:

`No map to state XX, you need to solve for more vectors in order to follow this state.`

This message usually means you need to increase the number of excited states to be calculated for a proper description of the one you’re interested in. Increase the number N for nstates=N in the route section at higher computational cost. A rule of thumb is to request at least 2 more states than the state of interest. This message can also reflect the fact that during the optimization the energy ordering changes between states, and can also mean that the ground state wave function is unstable, i.e., the energy of the excited state becomes lower than that of the ground state, in this case a single determinant approach is unviable and CAS should be used if the size of the molecule allows it. Excited state optimizations are tricky this way, in some cases the optimization may cross from one PES to another making it hard to know if the resulting geometry corresponds to the state of interest or another. Gaussian recommends changing the step size of the optimization from the default 0.3 Bohr radius to 0.1, but obviously this will make the calculation take longer.

`Opt=(MaxStep=10)`

If the minimum on the excited state potential energy surface (PES) doesn’t exist, then the excited state is not bound; take for example the first excited state of the H2 molecule which doesn’t show a minimum, and therefore the optimized geometry would correspond to both H atoms moving away from each other indefinitely (Figure 2). Nevertheless, a failed optimization doesn’t necessarily means the minimum does not exist and further analysis is required, for instance, checking the gradient is converging to zero while the forces do not. Figure 2. An unbound excited state with no minima ensures the dissociation of the system along the reaction coordinate

## Basis Set Superposition Error (BSSE). A short intro

Molecular Orbitals (MOs) are linear combinations of Atomic Orbitals (AOs), which in turn are linear combinations of other functions called ‘basis functions’. A basis, or more accurately a basis set, is a collection of functions which obey a set of rules (such as being orthogonal to each other and possibly being normalized) with which all AOs are constructed, and although these are centered on each atomic nucleus, the canonical way in which they are combined yield delocalized MOs; in other words, an MO can occupy a large space spanning several atoms at once. We don’t mind this expansion across a molecule, but what about between two molecules? Calculating the interaction energy between two or more molecular fragments leads to an artificial extra–stabilization term that stems from the fact that electrons in molecule 1 can occupy AO’s (or the basis functions which form them) centered on atoms from molecule 2.

Fundamentally, the interaction energy of any A—B dimer, Eint, is calculated as the energy difference between the dimer and the separately calculated energies for each component (Equation 1).

Eint = EAB – EA – EB (1)

However the calculation of Eint by this method is highly sensitive to the choice of basis set due to the Basis Set Superposition Error (BSSE) described in the first paragraph. The BSSE is particularly troublesome when small basis sets are used, due to the poor description of dispersion interactions but treating this error by just choosing a larger basis set is seldom useful for systems of considerable sizes. The Counterpoise method is a nifty correction to equation 1, in which EA and EB are calculated with the basis set of A and B respectively, i.e., only in EAB a larger basis set (that of A and B simultaneously) is used. The Counterpoise method calculates each component with the AB basis set (Equation 2)

EintCP = EABAB – EAAB– EBAB (2)

where the superscript AB means the whole basis set is used. This is accomplished by using ‘ghost‘ atoms with no nuclei and no electrons but empty basis set functions centered on them.

In Gaussian, BSSE is calculated with the Counterpoise method developed by Boys and Simon. It requires the keyword Counterpoise=N where N is the number of fragments to be considered (for an A—B system, N=2). Each atom in the coordinates list must be specified to which fragment pertains; additionally, the charge and multiplicity for each fragment and the whole supermolecular ensemble must be specified. Follow the example of this hydrogen fluoride dimer.

```%chk=HF2.chk
#P opt wB97XD/6-31G(d,p) Counterpoise=2

HF dimer

0,1 0,1 0,1
H(Fragment=1) 0.00 0.00 0.00
F(Fragment=1) 0.00 0.00 0.70
H(Fragment=2) 0.00 0.00 1.00
F(Fragment=2) 0.00 0.00 1.70```

For closed shell fragments the first line is straightforward but one must pay attention that the first pair of numbers in the charge multiplicity line correspond to the whole ensemble, whereas the folowing pairs correspond to each fragment in consecutive order. Fragments do not need to be specified contiguously, i.e., you don’t need to define all atoms for fragment 1 and after those the atoms for fragment 2, etc. They could be mixed and the program still assigns them correctly. Just as an example I typed wB97XD but any other method, DFT or ab initio, may be used; only semiempirical methods do not admit a BSSE calculation because they don’t make use of a basis set in the first place!

The output provides the corrected energy (in atomic units) for the whole system, as well as the BSSE correction (which added to the previous term yields the un-corrected energy of the system). Gaussian16 also provides these values in kcal/mol as ‘Complexation energies’ first raw (uncorrected) and then the corrected energy.

BSSE is always present and cannot be entirely eliminated because of the use of finite basis sets but it can be correctly dealt with if the Counterpoise method is included.

## Density Keyword in Excited State Calculations with Gaussian

I have written about extracting information from excited state calculations but an important consideration when analyzing the results is the proper use of the keyword density.

This keyword let’s Gaussian know which density is to be used in calculating some results. An important property to be calculated when dealing with excited states is the change in dipole moment between the ground state and any given state. The Transition Dipole Moment is an important quantity that allows us to predict whether any given electronic transition will be allowed or not. A change in the dipole moment (i.e. non-zero) of a molecule during an electronic transition helps us characterize said transition.

Say you perform a TD-DFT calculation without the density keyword, the default will provide results on the lowest excited state from all the requested states, which may or may not be the state of interest to the transition of interest; you may be interested in the dipole moment of all your excited states.

Three separate calculations would be required to calculate the change of dipole moment upon an electronic transition:

1) A regular DFT for the ground state as a reference
2) TD-DFT, to calculate the electronic transitions; request as many states as you need/want, analyze it and from there you can see which transition is the most important.
3) Request the density of the Nth state of interest to be recovered from the checkpoint file with the following route section:

`# TD(Read,Root=N) LOT Density=Current Guess=Read Geom=AllCheck`

replace N for the Nth state which caught your eye in step number 2) and LOT for the Level of Theory you’ve been using in the previous steps. That should give you the dipole moment for the structure of the Nth excited state and you can compare it with the one in the ground state calculated in 1). Again, if density=current is not used, only properties of N=1 will be printed.

## Orbital Contributions to Excited States

This is a guest post by our very own Gustavo “Gus” Mondragón whose work centers around the study of excited states chemistry of photosynthetic pigments.

When you’re calculating excited states (no matter the method you’re using, TD-DFT, CI-S(D), EOM-CCS(D)) the analysis of the orbital contributions to electronic transitions poses a challenge. In this post, I’m gonna guide you through the CI-singles excited states calculation and the analysis of the electronic transitions.

I’ll use adenine molecule for this post. After doing the corresponding geometry optimization by the method of your choice, you can do the excited states calculation. For this, I’ll use two methods: CI-Singles and TD-DFT.

The route section for the CI-Singles calculation looks as follows:

`%chk=adenine.chk%nprocshared=8%mem=1Gb#p CIS(NStates=10,singlets)/6-31G(d,p) geom=check guess=read scrf=(cpcm,solvent=water)adenine excited states with CI-Singles method0 1--blank line--`

I use the same geometry from the optimization step, and I request only for 10 singlet excited states. The CPCP implicit solvation model (solvent=water) is requested. If you want to do TD-DFT, the route section should look as follows:

`%chk=adenine.chk%nprocshared=8%mem=1Gb#p FUNCTIONAL/6-31G(d,p) TD(NStates=10,singlets) geom=check guess=read scrf=(cpcm,solvent=water)adenine excited states with CI-Singles method0 1--blank line--`

Where FUNCTIONAL is the DFT exchange-correlation functional of your choice. Here I strictly not recommend using B3LYP, but CAM-B3LYP is a noble choice to start.

Both calculations give to us the excited states information: excitation energy, oscillator strength (as f value), excitation wavelength and multiplicity:

Excitation energies and oscillator strengths:

` Excited State   1:      Singlet-A      6.3258 eV  196.00 nm  f=0.4830  <S**2>=0.000      11 -> 39        -0.00130      11 -> 42        -0.00129      11 -> 43         0.00104      11 -> 44        -0.00256      11 -> 48         0.00129      11 -> 49         0.00307      11 -> 52        -0.00181      11 -> 53         0.00100      11 -> 57        -0.00167      11 -> 59         0.00152      11 -> 65         0.00177`

The data below corresponds to all the electron transitions involved in this excited state. I have to cut all the electron transitions because there are a lot of them for all excited states. If you have done excited states calculations before, you realize that the HOMO-LUMO transition is always an important one, but not the only one to be considered. Here is when we calculate the Natural Transition Orbitals (NTO), by these orbitals we can analyze the electron transitions.

For the example, I’ll show you first the HOMO-LUMO transition in the first excited state of adenine. It appears in the long list as follows:

35 -> 36         0.65024

The 0.65024 value corresponds to the transition amplitude, but it doesn’t mean anything for excited state analysis. We must calculate the NTOs of an excited state from a new Gaussian input file, requesting from the checkpoint file we used to calculate excited states. The file looks as follows:

`%Oldchk=adenine.chk%chk=adNTO1.chk%nproc=8%mem=1Gb#p SP geom=allcheck guess=(read,only) density=(Check,Transition=1) pop=(minimal,NTO,SaveNTO)`

I want to say some important things right here for this last file. See that no level of theory is needed, all the calculation data is requested from the checkpoint file “adenine.chk”, and saved into the new checkpoint file “adNTO1.chk”, we must use the previous calculated density and specify the transition of interest, it means the excited state we want to analyze. As we don’t need to specify charge, multiplicity or even the comment line, this file finishes really fast.

After doing this last calculation, we use the new checkpoint file “adNTO1.chk” and we format it:

`formchk -3 adNTO1.chk adNTO1.fchk`

If we open this formatted checkpoint file with GaussView, chemcraft or the visualizer you want, we will see something interesting by watching he MOs diagram, as follows:

We can realize that frontier orbitals shows the same value of 0.88135, which means the real transition contribution to the first excited state. As these orbitals are contributing the most, we can plot them by using the cubegen routine:

`cubegen 0 mo=homo adNTO1.fchk adHOMO.cub 0 h`

This last command line is for plotting the equivalent as the HOMO orbital. If we want to plot he LUMO, just change the “homo” keyword for “lumo”, it doesn’t matter if it is written with capital letters or not.

You must realize that the Natural Transition Orbitals are quite different from Molecular Orbitals. For visual comparisson, I’ve printed also the molecular orbitals, given from the optimization and from excited states calculations, without calculating NTOs:

These are the molecular frontier orbitals, plotted with Chimera with 0.02 as the isovalue for both phase spaces:

The frontier NTOs look qualitatively the same, but that’s not necessarily always the case:

If we analyze these NTOs on a hole-electron model, the HOMO refers to the hole space and the LUMO refers to the electron space.

Maybe both orbitals look the same, but both frontier orbitals are quite different between them, and these last orbitals are the ones implied on first excited state of adenine. The electron transition will be reported as follows:

If I can do a graphic summary for this topic, it will be the next one:

NTOs analysis is useful no matter if you calculate excited states by using CIS(D), EOM-CCS(D), TD-DFT, CASSCF, or any of the excited states method of your election. These NTOs are useful for population analysis in excited states, but these calculations require another software, MultiWFN is an open-source code that allows you to do this analysis, and another one is called TheoDORE, which we’ll cover in a later post.

## Natural Transition Orbitals (NTOs) Gaussian

The canonical molecular orbital depiction of an electronic transition is often a messy business in terms of a ‘chemical‘ interpretation of ‘which electrons‘ go from ‘which occupied orbitals‘ to ‘which virtual orbitals‘.

Natural Transition Orbitals provide a more intuitive picture of the orbitals, whether mixed or not, involved in any hole-particle excitation. This transformation is particularly useful when working with the excited states of molecules with extensively delocalized chromophores or multiple chromophoric sites. The elegance of the NTO method relies on its simplicity: separate unitary transformations are performed on the occupied and on the virtual set of orbitals in order to get a localized picture of the transition density matrix.

 R. L. Martin, J. Chem. Phys., 2003, DOI:10.1063/1.1558471.

In Gaussian09:
After running a TD-DFT calculation with the keyword TD(Nstates=n) (where n = number of states to be requested) we need to take that result and launch a new calculation for the NTOs but lets take it one step at a time. As an example here’s phenylalanine which was already optimized to a minimum at the B3LYP/6-31G(d,p) level of theory. If we take that geometry and launch a new calculation with the TD(Nstates=40) in the route section we obtain the UV-Vis spectra and the output looks like this (only the first three states are shown):

```Excitation energies and oscillator strengths:

Excited State 1: Singlet-A 5.3875 eV 230.13 nm f=0.0015 <S**2>=0.000
42 -> 46 0.17123
42 -> 47 0.12277
43 -> 46 -0.40383
44 -> 45 0.50838
44 -> 47 0.11008
This state for optimization and/or second-order correction.
Total Energy, E(TD-HF/TD-KS) = -554.614073682
Copying the excited state density for this state as the 1-particle RhoCI density.

Excited State 2: Singlet-A 5.5137 eV 224.86 nm f=0.0138 <S**2>=0.000
41 -> 45 -0.20800
41 -> 47 0.24015
42 -> 45 0.32656
42 -> 46 0.10906
42 -> 47 -0.24401
43 -> 45 0.20598
43 -> 47 -0.14839
44 -> 45 -0.15344
44 -> 47 0.34182

Excited State 3: Singlet-A 5.9254 eV 209.24 nm f=0.0042 <S**2>=0.000
41 -> 45 0.11844
41 -> 47 -0.12539
42 -> 45 -0.10401
42 -> 47 0.16068
43 -> 45 -0.27532
43 -> 46 -0.11640
43 -> 47 0.16780
44 -> 45 -0.18555
44 -> 46 -0.29184
44 -> 47 0.43124```

The oscillator strength is listed on each Excited State as “f” and it is a measure of the probability of that excitation to occur. If we look at the third one for this phenylalanine we see f=0.0042, a very low probability, but aside from that the following list shows what orbital transitions compose that excitation and with what energy, so the first line indicates a transition from orbital 41 (HOMO-3) to orbital 45 (LUMO); there are 10 such transitions composing that excitation, visualizing them all with canonical orbitals is not an intuitive picture, so lets try the NTO approach, we’re going to take excitation #10 for phenylalanine as an example just because it has a higher oscillation strength:

```%chk=Excited State 10: Singlet-A 7.1048 eV 174.51 nm f=0.3651 <S**2>=0.000
41 -> 45 0.35347
41 -> 47 0.34685
42 -> 45 0.10215
42 -> 46 0.17248
42 -> 47 0.13523
43 -> 45 -0.26596
43 -> 47 -0.22995
44 -> 46 0.23277```

Each set of NTOs for each transition must be calculated separately. First, copy you filename.chk file from the TD-DFT result to a new one and name it after the Nth state of interest as shown below (state 10 in this case). NOTE: In the route section, replace N with the number of the excitation of interest according to the results in filename.log. Run separately for each transition your interested in:

```#chk=state10.chk

#p B3LYP/6-31G(d,p) Geom=AllCheck Guess=(Read,Only) Density=(Check,Transition=N) Pop=(Minimal,NTO,SaveNTO)

0 1
--blank line--```

By requesting SaveNTO, the canonical orbitals in the state10.chk file are replaced with the NTOs for the 10th excitation, this makes it easier to plot since most visualizers just plot whatever set of orbitals they read in the chk file but if they find the canonical MOs then one would need to do some re-processing of them. This is much more straightforward.

Now we format our chk files into fchk with the formchk utility:

`formchk -3 filename.chk filename.fchkformchk -3 state10.chk state10.fchk`

If we open filename.fchk (the file where the original TD-DFT calculation is located) with GaussView we can plot all orbitals involved in excited state number ten, those would be seven orbitals from 41 (HOMO-3) to 47 (LUMO+2) as shown in figure 1. Figure 1. Canonical orbitals involved in the 10th excited state according to the TD-DFT calculation

If we now open state10.fchk we see that the numbers at the side of the orbitals are not their energy but their occupation number particular to this state of interest, so we only need to plot those with highest occupations, in our example those are orbitals 44 and 45 (HOMO and LUMO) which have occupations = 0.81186; you may include 43 and 46 (HOMO-1 and LUMO+1, respectively) for a much more complete description (occupations = 0.18223) but we’re still dealing with 4 orbitals instead of 7. Figure 2. Natural Transition Orbitals for Phenylalanine. Orbital 44 (particle) and Orbital 45 (hole) exhibit the largest occupations for Excited State No. 10

The NTO transition 44 -> 45 is far easier to conceptualize than all the 10 combinations given in the canonical basis from the direct TD-DFT calculation. TD-DFT provides us with the correct transitions, NTOs just paint us a picture more readily available to the chemist mindset.

NOTE: for G09 revC and above, the %OldChk option is available, I haven’t personally tried it but using it to specify where the excitations are located and then write the NTOs of interest into a new chk file in the following way, thus eliminating the need of copying the original chk file for each state:

`%OldChk=filename.chk%chk=stateN.chk`

NTOs are based on the Natural Hybrid orbitals vision by Löwdin and others, and it is said to be so straightforward that it has been re-discovered from time to time. Be that as it may, the NTO visualization provides a much clearer vision of the excitations occurring during a TD calculation.

Thanks for reading, stay home and stay safe during these harsh days everyone. Please share, rate and comment this and other posts.

## Useful Thermochemistry from Gaussian Calculations

Statistical Mechanics is the bridge between microscopic calculations and thermodynamics of a particle ensemble. By means of calculating a partition function divided in electronic, rotational, translational and vibrational functions, one can calculate all thermodynamic functions required to fully characterize a chemical reaction. From these functions, the vibrational contribution, together with the electronic contribution, is the key element to getting thermodynamic functions.

Calculating the Free Energy change of any given reaction is a useful approach to asses their thermodynamic feasibility. A large negative change in Free Energy when going from reagents to products makes up for a quantitative spontaneous (and exothermic) reaction, nevertheless the rate of the reaction is a different story, one that can be calculated as well.

Using the freq option in your route section for a Gaussian calculation is mandatory to ascertain the current wave function corresponds to a minimum on a potential energy hypersurface, but also yields the thermochemistry and thermodynamic values for the current structure. However, thermochemistry calculations are not restricted to minima but it can also be applied to transition states, therefore yielding a full thermodynamic characterization of a reaction mechanism.

A regular freq calculation yields the following output (all values in atomic units):

```Zero-point correction=                           0.176113 (Hartree/Particle)
Thermal correction to Energy=                    0.193290
Thermal correction to Enthalpy=                  0.194235
Thermal correction to Gibbs Free Energy=         0.125894
Sum of electronic and zero-point Energies=           -750.901777
Sum of electronic and thermal Energies=              -750.884600
Sum of electronic and thermal Enthalpies=            -750.883656
Sum of electronic and thermal Free Energies=         -750.951996```

For any given reaction say A+B -> C one could take the values from the last row (lets call it G) for all three components of the reaction and perform the arithmetic: DG = GC – [GA + GB], so products minus reagents.

By default, Gaussian calculates these values (from the previously mentioned partition function) using normal conditions, T = 298.15 K and P = 1 atm. For an assessment of the thermochemistry at other conditions you can include in your route section the corresponding keywords Temperature=x.x and Pressure=x.x, in Kelvin and atmospheres, respectively.

(Huge) Disclaimer: Although calculating the thermochemistry of any reaction by means of DFT calculations is a good (and potentially very useful) guide to chemical reactivity, getting quantitative results require of high accuracy methods like G3 or G4 methods, collectively known as Gn mehtods, which are composed of pre-defined stepwise calculations. The sequence of these calculations is carried out automatically; no basis set should be specified. Other high accuracy methods like CBS-QB3 or W1U can also be considered whenever Gn methods are too costly.

## Atom specifications unexpectedly found in input stream.

“Well, where else were they supposed to appear?”

I was sent this error along with the previous question for a failed optimization. Apparently there is no answer in the internet (I quickly checked) so here it is:

Gaussian is confused about finding atomic coordinates because there is also a geom=check instruction placed in the route section, i.e., it was told to retrieve the atomic coordinates from a checkpoint and then it was given those atomic coordinates within the input so it doesn’t know what you mean and exits.