Category Archives: Visualization
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.
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.
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.fchk
formchk -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.
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.
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:
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.
We’ve covered some common errors when dealing with formatted checkpoint files (*.fchk) generated from Gaussian, specially when analyzed with the associated GaussView program. (see here and here for previous posts on the matter.)
Prof. Neal Zondlo from the University of Delaware kindly shared this solution with us when the following message shows up:
CConnectionGFCHK::Parse_GFCHK() Missing or bad data: Rbond Line Number 1234
The Rbond label has to do with the connectivity displayed by the visualizer and can be overridden by close examination of the input file. In the example provided by Prof. Zondlo he found the following line in the connectivity matrix of the input file:
2 9 0.0
which indicates a zero bond order between atoms 2 and 9, possibly due to their proximity. He changed the line to simply
So editing the connectivity of your atoms in the input can help preventing the Rbond message.
I hope this helps someone else.
A couple of weeks ago I posted a solution for a common error regarding .fchk files that will display the error below when opened with GaussView5.0. As I expected, this error has to do with the use of diffuse functions in the basis set and is related to a change of format between Gaussian versions.
CConnectionGFCHK::Parse_GFCHK() Missing or bad data: Alpha Orbital Energies Line Number 1234
Although the method described in the previous post works just fine, the following update is a better approach. Due to a change of spelling between G03 and G09 (which has been corrected for G09 but not available for GV versions prior to 5.0.9) one must change “independent” for “independant”
To make the change directly from the terminal the following command is needed:
sed -i 's/independent/independant/g' file.fchk
Alternatively you can redirect the output to a new file
sed -e 's/independent/independant/g' file.fchk > newfile.fchk
if you want to keep the old version and work with a new one.
Of course this edition can be performed manually with any text editor available (for example if you work in Windows) but solutions from the terminal always seem easier and a lot more fun to me.
Thanks to Dr. Fernando Cortés for sharing his insight into this issue.
I’ve found the following error regarding the opening of .fchk files in GaussView5.0.
CConnectionGFCHK::Parse_GFCHK() Missing or bad data: Alpha Orbital Energies Line Number 1234
The error is prevented to a first approximation (i.e. it at least will allow GV to open and visualize the file but other issues may arise) by opening the file and modifying the number of basis functions to equal the number of independent functions (which is lower)
FILE HEADER FOpt RM062X 6-311++G(d,p) Number of atoms I 75 Info1-9 I N= 9 163 163 0 0 0 110 2 18 -502 Charge I 0 Multiplicity I 1 Number of electrons I 314 Number of alpha electrons I 157 Number of beta electrons I 157 Number of basis functions I 1199 Number of independent functions I 1199 Number of point charges in /Mol/ I 0 Number of translation vectors I 0 Atomic numbers I N= 75 ... ... ... ...
Once both numbers match you can open the file normally and work with it. My guess is this will continue to happen with highly polarized basis sets but I need to run some tests.
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.
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.
> It’s possible to open a formated checkpoint file, containing the NBOs, in chemcraft.
%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!
This is the first time I reblog a post from a fellow computational chemist and the reason why I do it is because of its beautiful simplicity and usefulness. Given the scope this blog has taken I think this post becomes most appropriate. This post will show you how to create an energy level diagram using nothing but MS Excel.
Kudos to ‘Eutactic’, from Australia, for coming up with a nice solution to this problem. Check out his blog at eutactic.wordpress.com.
Thanks for letting me repost it 🙂
I worked out a very quick and easy way to generate level schemes in Excel, based on a query from one of the other students in the group. Normally I would resort to something like the astonishing TikZ for this sort of task, however our group is very much a Microsoft Office ‘What You See Is A Metaphor For Cosmic Horror‘ group and recommending that a colleague learns two new markup languages to produce a figure is probably not helpful in the short term. One of the issues with charting energy levels in Excel is that levels are typically represented by horizontal bars connected at their vertices with lines representing transitions. Whilst Excel does have a horizontal bar as a marker, it possesses two show-stopping limitations:
- It is only uniformly scalable, and can only be scaled so far – we cannot make it anywhere near wide and…
View original post 222 more words
Due to extensive popular demand, I hereby make available the necessary files to run Molekel in its old 4.3 version. The program has been compiled to work under Windows 32 bit architecture. Just extract it and place the main folder (as provided here) in any location and run the .exe file located inside. You can generate a direct access to it from your desktop and it even includes a small icon to be used for this purpose.
Also, the manual in pdf format is included; in it you can find the proper citation which must be included in any publication that makes use of Molekel. Just in case you can’t find it, here it is:
MOLEKEL 4.3, P. Flükiger, H.P. Lüthi, S. Portmann, J. Weber, Swiss Center for Scientific
Computing, Manno (Switzerland), 2000-2002.
Stefan Portmann & Hans Peter Lüthi. MOLEKEL: An Interactive Molecular Graphics Tool.
CHIMIA (2000) 54 766-770.
Now some considerations:
- This is an old program. It was generated back in the WindowsXP days, so it wouldn’t be a surprise if it doesn’t work in more recent platforms or under any other 32 bit OS.
- The manual is included. Please read it. This blog is not Molekel’s help desk; I may help but I can’t solve everything, I just don’t have the time for it.
- I didn’t participate/collaborate/helped or got involved in any way in the development of this program, i.e., don’t shoot the messenger. The Molekel homepage is: http://molekel.cscs.ch/wiki/pmwiki.php
- I strongly recommend to make use of the NEW version. A bit more obscure but also great once you figure it out, plus there is available support for it from the actual developers.
- I also strongly recommend to look over the internet for other visualization softwares. I don’t recall having reviewed any in this blog. Perhaps some other time.
- Since this is not my development I will remove it from the server upon the request of the rightful owners A.S.A.P! My guess is they wont mind all that much since its an old version and it was given away for free from their server anyway.
- I can’t think of anything else to put on this list right now but I reserve the right to come back to it and add something more. I just don’t want any trouble.
So, here it is! Right click on the link and download it; Use it to generate nice plots of your orbitals, densities, electrostatic potentials, etc. Consider this a Happy New Year’s gift!
Rate and comment this and all the other posts you find interesting in this blog. Please!
UPDATE: Thanks to Yuekui Wang for the following information.
This copy doen’t work on some WinXP machine with ATI monitor card. The original copy is still available on the cscs web site. Download link is as follows:
It works fine on many macjines, I am sure.