# Category Archives: NBO

## 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.

## 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.

## Python scripts for calculating Fukui Indexes

One of the most popular posts in this blog has to do with calculating Fukui indexes, however, when dealing with a large number of molecules, our described methodology can become cumbersome since it requires to manually extract the population analysis from two or three different output files and then performing the arithmetic on them separately with a spreadsheet or something.

Our new team member Ricardo Loaiza has written a python script that takes the three aforementioned files and yields a .csv file with the calculated Fukui indexes, and it even points out which of the atoms exhibit the largest values so if you have a large molecule you don’t have to manually check for them. We have also a batch version which takes all the files in any given directory and performs the Fukui calculations for each, provided it can find file triads with the naming requirements described below.

Output files must be named filename.log (the N electrons reference state), filename_plus.log (the state with N+1 electrons) and filename_minus.log (the N-1 electrons state). Another restriction is that so far these scripts only work with NBO population analysis as provided by the NBO3.1 program available in the various versions of Gaussian. I imagine the listing is similar in NBO5.x and NBO6.x and so it should work if you do the population analysis with them.

The syntax for the single molecule version is:

`python fukui.py filename.log filename_minus.log filename_plus.log`

For the batch version is:

`./fukuiPorLote.sh`

(Por Lote means In Batch in Spanish.)

These scripts are available via GitHub. We hope you find them useful, and you do please let us know whether here at the comments section or at our GitHub site.

## Natural Bond Orbitals (NBO) Visualization with Chemcraft

It’s been a long time since I last posted something and so many things have happened in our research group! I should catch up with them in short but times have just been quite hectic.

I’m glad to publicly thank Prof. Frank Weinhold’s gesture to include this blog in the bibliography section of the new NBO6.0 website under the NBO-Related Websites tab.

Here is a contribution from Igor Marques at the University of Aveiro in Portugal (Group Website); the original text can be found as a comment in the original NBO Visualization post but it is pretty much the same thing you can find in this post. Here is a link to Chemcraft’s website. Thanks for sharing this, Igor!

=> Examples provided by Igor Marques used Chemcraft Version 1.7, build 365 <=

In the Gaussian input, with the NBORead option included under the population keyword, we should include the PLOT option as illustrated below. The gfoldprint keyword will print the basis set to the output file in the old G03 format. Some visualization programs require a certain format of the basis set to be printed to the output file in order to plot orbitals and other surfaces like the electron density; therefore, if you want to play safe, use gfoldprint, gfprint and gfinput in the same line. gfprint will print the basis set as a list but in the new G09 format, whereas gfinput will print the basis set using Gaussian’s own input format. (The used level of theory and number of shared processors are shown as illustrations only; also the Opt keyword is not fundamental to the visualization of the NBO’s)

```%chk=filename.chk
%nprocshared=8

#P b3lyp/6-311++g** Opt pop=(full,nboread) gfoldprint
filename

0 1
molecular coordinates
\$NBO BNDIDX PLOT \$END```

this will generate files from *.31 to *.41
For the visualization of NBOs, you’ll need FILE.31 and FILE.37. Open FILE.31 from chemcraft. It will automatically detect FILE.37 (if in the same directory).

`Tools > Orbitals > Render molecular orbitals`

select the NBOs of interest (whcih are in the same order of the output),

`Adjust settings > OK`

On the left side of the window, select the NBO of interest and then click on ‘show isosurface’. Adjust the remaining settings. To represent another orbital, click on ‘keep this surface’ and then select another orbital from the rendered set and follow the previous steps.

Some Considerations:

> It’s possible to open a formated checkpoint file, containing the NBOs, in chemcraft.
Gaussian input:

```%Chk=filename.chk
%nprocshared=4
#P b3lyp/6-311++g** Opt pop=(full,nboread,savenbo) gfoldprint
filename

0 1
molecular coordinates
\$NBO BNDIDX \$END```

the procedure is identical, but it is only necessary to read the *fchk file and then render the desired orbitals.
However, two problems might arise:
a) Orbitals in the checkpoint are reordered, thus requiring some careful inspection of the output.
b) Sometimes, for a larger molecule, the checkpoint might not be properly saved and the Gaussian job (as previously reported – http://goo.gl/DrSgA ) will end with:

Failed in SchOr1 in NBStor.
Error termination via Lnk1e in /data/programs/g09/l607.exe at Wed Mar 6 15:27:33 2013.

****

As usual, thanks to all for reading/commenting/rating this and other posts in this blog!