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:
#p CIS(NStates=10,singlets)/6-31G(d,p) geom=check guess=read scrf=(cpcm,solvent=water)
adenine excited states with CI-Singles method
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:
#p FUNCTIONAL/6-31G(d,p) TD(NStates=10,singlets) geom=check guess=read scrf=(cpcm,solvent=water)
adenine excited states with CI-Singles method
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:
#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.
It seems a bit weird that there isn’t much information on this topic on the internet. Recently I’ve had to calculate some of these indices to explain an anomalous behaviour in lactones formation and out of curiosity I ran a small search on the net about how to calculate them. Most of the information I retrieved were papers dealing with calculated Fukui and condensed Fukui indices, but unless you count with electronic subscriptions to the corresponding journals you were left in the dark. Moreover, even if you get to read the paper they will tell you as much about how to calculate Fukui indices as they tell you about the Hartree-Fock procedure details.
Therefore here I post this information specifically under “how to calculate Fukui indices” so others might find it. It seems to me this blog is taking a rather educational turn which was not its original intention. Still, if I can atract people to read about my work while finding useful information I’m glad to do it.
Fukui indices are, in short, reactivity indices; they give us information about which atoms in a molecule have a larger tendency to either loose or accept an electron, which we chemist interpret as which are more prone to undergo a nucleophilic or an electrophilic attack, respectively. This in turn has to do with a molecule’s tendency of becoming polarized in the presence of an external field or upon the change of electron density. The key word here is electron density, the whole idea behind Fukui’s lies in the realm of conceptual DFT.
Fukui functions are defined as the functional derivative of the chemical potential respect to the external potential (the one produced by the nuclei) at a constant electron number. Since the chemical potential is defined as the derivative of the density functional respect to the electron density, fukui functions are also defined as the derivative of the electron density respect to the number of electrons at a constant potential, and this latter definition is what we want to work with because it means that we can calculate how the density changes at every point (since it is already different at every point r) when adding or removing an electron while keeping the potential constant (that is the position of the nuclei, in other words mantaining the molecular geometry). The name fukui function stems from the fact that these added/removed electrons go into the frontier or Fukui orbitals HOMO/LUMO, but the reality is that the definition was conceived by Yang and Parr (like almost everything in conceptual DFT).
To practically perform these calculations the finite differences method is employed and so the condensed to atom fukui index is obtained.
Electrophilicity of atom A in molecule M (of N electrons)
fA+ = PA(N+1)-PA(N)
Nucleophilicity of atom A in molecule M (of N electrons)
fA– = PA(N)-PA(N-1)
Radical attack susceptibility of atom A in molecule M (of N electrons)
fA0 = 1/2[PA(N+1)-PA(N-1)]
where P stands for the population of atom A in molecule M. If you want to analyze the fukui function at an ionized species the procedure is the same but most people need to beware that N is the number of electrons of the original ion! (i.e. the species you are trying to analyze), sometimes it may be confusing if you are trying to analyze the nucleophilicity (f-) of an atom in a cationic species (M+).
The population analysis on the + or – system has to be performed at the same equilibrium geometry as the original molecule! If we optimize again the system then we are letting the system relax and therefore we loose information on the polarization of the electron density upon the change in number of electrons.
“Negative indices are meaningless and should be disregarded” This is a common statement that ever since the definition of Fukui indices has been regarded as true; however there have been a couple of recent publications (see below) that defy to some extent such notion.
Major drawback: Since we are ultimately dealing with occupation numbers on each atom these indices are very sensitive towards changes in basis sets and population analysis paradigm. I strongly recomend never to take those numbers as absolutes, only as comparative parameters within the same system! I also find useful to compute them using more than just one method and one level of theory in order to further confirm or even to dismiss the trends observed. Natural population analysis and AIM are much more robust than simple Mulliken PA.
Hopefully this will be helpful to people trying to calculate fukui functions and also to understand what they mean. Subscribe to this post for further updates (you never know). Rates and comments are always most welcome specially if you found this post interesting or useful. Cheers!
Further reading.-Computational Chemistry by Errol Lewars
-Bultinck et al. Negative Fukui functions: New insights based on electronegativity
equalization J.Chem.Phys 118 (10) 4349-4356 (2003)
-Melin J, Ayers PW, Ortiz JV. Removing electrons can increase the electron density: a computational study of negative Fukui functions. J Phys Chem A. 2007 111(40):10017-9
-Bultnick & Cabo-Dorca Negative and Infinite Fukui Functions: The Role of Diagonal Dominance in the Hardness Matrix J. Mat. Chem. 34 (2003) 67-74
PLEASE CHECK THIS OTHER POST FOR A NICE SCRIPT OF OURS TO AUTOMATE THIS PROCEDURE. THANKS!