looking at molecules 1

this is the first part of an upcoming multi-post guide, in which I want to describe how to address photochemical questions about a certain molecule or class of molecules. Therefore, I think it is a good approach to first line out the process of the quantum chemical investigation in a step-like structure and then go into details of each step in future posts. So lets start from the beginning:

For starters, I want to get an overview about the molecules inner (electronic) structure and reactivity. But before I turn on the electric heating of the server room I thoroughly consult existing literature and any available scientists about reactivity and photochemical behaviour of the molecule(s). This helps to avoid getting lost in the vast amount of information that the quantum-chemical machinery can provide you with.

Secondly, on the basis of the literature and my experience I hypothesize about the electronic structure of the molecule and only afterwards start my calculations. At this point, you inevitably face the choice of the methodology. Depending on the size of the system, I usually start off with some solid double-zeta (SVP) ab-initio treatment (if doable CCSD, otherwise SCS/MOS-MP2) and double-hybrid functional (B2-PLYP). At this point, the calculations are limited to ground state geometry of the molecule and I try to get familiar with the MOs. To assure reliability, it is always a good idea to compare the outcome various methods and learn from the differences in the predictions. This way you can  e.g. identify highly correlated states whose energy changes a lot when going from second- to third-order methods or charge-transfer states when using TDDFT.


The third part of my guide will be concerned with the identification of the relevant coordinates of the system under investigation. Since most photochemical processes involve multi-state degeneracies (e.g. conical intersections and/or intersystem crossing), a key step in any quantum-chemical investigations is to learn along which geometrical parameters the involves states exhibit a different run and intersect. This can be intra-molecular coordinates such as bond length and angles or an inter-molecular distance coordinate along which a charge-transfer state of a molecular complex intersects with locally excited state.

The fourth and last section is about getting the numbers right. While in the first three parts we were satisfied with an approximate treatment of our system, we now want to get as quantitative as possible. Therefore, you should have some experimental numbers at hand to compare-to and a thorough insight into the electronic structure in the system as well as the experiments.



Now that I laid out the tasks, the next posts will be dealing with the first of those steps.

Cheers,
Jan

On the molecular mechanism of non-radiative decay in nitrobenzene and the unforseen molecular challenges this small molecule holds for electronic structure theory

This lengthy title is from our recently published paper that had been mentioned earlier in this blog. Amongst with some other things, this paper or rather the work that went into it also was the reason for the long pause since the last post. Now that it finally completed the process of publication, I will continue with a tutorial on how I look at molecules with the tools of a theoretical chemist.
There will be multiple posts on this topic in the upcoming weeks, but first, here is the link:

http://pubs.rsc.org/en/content/articlelanding/2014/cp/c4cp01232a

So long!

Shiny Pictures

Today, Felix introduced me to his new density analysis code and how he translates the resulting .cube-files into shiny 3D-pictures. 

His density analysis stuff comes in really handy for the characterization and interpretation of excited-states. There will be a paper out very soon, which you should definitely read (and cite)!

Transition densities of the bright
2 A1 (top), 1 B1 (middle), and
lowest 1 A2 singlet excited states
of nitrobenzene at the ADC2/
cc-pVDZ level of theory

What impressed me even more coming in even handier was the degree of automation he established for turning the .cube-files into pictures via VMD. He wrote a post about that some time ago and since then I guess I could have saved 2 or 3 days of work if I read and used it ever since. If you belong to the group of people using VMD to render MOs the manual way, you should really look into this his post or continue reading here.

After he showed me his work today, I basically stole his bash script, added options to render either two or three isosurfaces instead of just one and then spent the whole day adjusting the surface properties to make the pictures as shiny and informative as possible. If you like the results shown on the right side (using 3 isosurfaces with the standard isovalues of our script (0.0128, 0.0032, 0.0008), you might want to read on and learn how the script works:

If executed in a folder, it iterates over all .cube files present to write out 4 files:
load_all_plt.vmd
plot_all.vmd
convert.bash
vmd_plots.html.

Next, the user reads the geometry (e.g. from any xyz file) into vmd (vmd file.xyz), which needs to have the same orientation as in the job that created the cube files.

Afterwards, the vmd option "load visualization state" is used to read in the "load_all_plt.vmd" file. This script basically tells VMD to load all the .cube files and sets the options for the surfaces. At this point, the user can adapt the visualization and orientation of the molecule according to his wishes.

After selecting a reasonable orientation for the molecule in the usual way (with the mouse), the "plot_all.vmd" file is loaded, causing VMD to successively render all .cube files found in the folder yielding a bunch of .tga picture files.

Now, the convert.bash script is executed from the terminal to turn all those .tga files into .png files, which can be viewed in a very convenient way - by just opening the .html file with any browser. - DONE

Again, many thanks to Felix for this script! Thank him by (using his analysis tools and) citing is work!

And here is the script:
#!/bin/bash
# 0. $pointval mo 65-74 (Turbomole), plots (QChem), ...
# 1. call this script
# 2. open the molecular structure file in VMD
# 3. load the .plt/.cube files and some settings
#    - "Load state" load_all_plt.vmd
#    - click "Apply" in "Graphical Representations"
# 4. adjust perspective
# 5. "Load state" plot_all.vmd

###
ifmt=cube
ofmt=tga
out=load_all_plt.vmd
plot=plot_all.vmd
conv=convert.bash
html=vmd_plots.html
ncol=4
###

echo 'USAGE: $0 [<2 or 3 surfaces, STD = 3 >] [<highest iso, STD(2/3) = (0.1/0.128)>]'

isotemp=$2

if [ $1 -eq 2 ]
then
  isov=${isotemp:=0.01}
  isov2=`echo $isov'/8' | bc -l`
  isov3=0.99

#isovalue of 0.99 will produce no surface
  echo 'Using 2 surfaces for isovalues:'
  echo $isov $isov2
  echo "material change opacity Glass3 0.150000" > $out
  echo "material change diffuse Glass3 0.10000" >> $out
elif [ $1 -eq 3 ]
then
  isov=${isotemp:=0.0128}
  isov2=`echo $isov'/4' | bc -l`
  isov3=`echo $isov'/16' | bc -l`
  echo 'Using 3 surfaces for isovalues:'
  echo $isov $isov2 $isov3
  echo "material change opacity Glass3 0.400000" > $out
else
   echo "Please enter 2 or 3 for # of surfaces!"
   echo "Falling back to 3 surfaces with standard values"
   isov=0.0128
   isov2=0.0032
   isov3=0.0008
  echo "Using standard isovalues: "
  echo "Values: "$isov $isov2 $isov3
  echo "material change opacity Glass3 0.400000" > $out
fi

echo "axes location Off" >> $out
echo "display projection Orthographic" >> $out
echo "display rendermode GLSL" >> $out
echo "display depthcue off" >> $out
echo "color Display Background white" >> $out
echo "menu graphics on" >> $out
echo "material change diffuse Ghost 0.000000" >> $out
echo "material change ambient Ghost 0.300000" >> $out
echo "material change opacity Ghost 0.100000" >> $out
echo "material change shininess Ghost 0.000000" >> $out
echo "mol addrep 0" >> $out
echo "mol addrep 0" >> $out
echo "mol addrep 0" >> $out
echo "mol addrep 0" >> $out
echo "mol addrep 0" >> $out
echo "mol addrep 0" >> $out
echo "mol modmaterial 1 0 Opaque" >> $out
echo "mol modmaterial 2 0 Opaque" >> $out
echo "mol modmaterial 3 0 Glass3" >> $out
echo "mol modmaterial 4 0 Glass3" >> $out
echo "mol modmaterial 5 0 Ghost" >> $out
echo "mol modmaterial 6 0 Ghost" >> $out
echo "mol modstyle 1 0 Isosurface  $isov 0 0 0 1 1" >> $out
echo "mol modstyle 2 0 Isosurface -$isov 0 0 0 1 1" >> $out
echo "mol modstyle 3 0 Isosurface  $isov2 0 0 0 1 1" >> $out
echo "mol modstyle 4 0 Isosurface -$isov2 0 0 0 1 1" >> $out
echo "mol modstyle 5 0 Isosurface  $isov3 0 0 0 1 1" >> $out
echo "mol modstyle 6 0 Isosurface -$isov3 0 0 0 1 1" >> $out
echo "mol modcolor 1 0 ColorID 0" >> $out
echo "mol modcolor 2 0 ColorID 1" >> $out
echo "mol modcolor 3 0 ColorID 0" >> $out
echo "mol modcolor 4 0 ColorID 1" >> $out
echo "mol modcolor 5 0 ColorID 0" >> $out
echo "mol modcolor 6 0 ColorID 1" >> $out
echo "" > $plot

echo "#!/bin/bash" > $conv
chmod +x $conv

echo -e "<html>\n<head></head>\n<body>" > $html
echo -e "<table>\n<tr>" >> $html

N=0
for I in *$ifmt
do
   echo "mol addfile $I" >> $out
   echo "mol modstyle 1 0 Isosurface  $isov $N 0 0 1 1" >> $plot
   echo "mol modstyle 2 0 Isosurface -$isov $N 0 0 1 1" >> $plot
   echo "mol modstyle 3 0 Isosurface  $isov2 $N 0 0 1 1" >> $plot
   echo "mol modstyle 4 0 Isosurface -$isov2 $N 0 0 1 1" >> $plot
   echo "mol modstyle 5 0 Isosurface  $isov3 $N 0 0 1 1" >> $plot
   echo "mol modstyle 6 0 Isosurface -$isov3 $N 0 0 1 1" >> $plot
   echo "render TachyonInternal $I.$ofmt" >> $plot

   echo "convert $I.$ofmt $I.png" >> $conv
   echo "rm $I.$ofmt" >> $conv

   echo "<td><img src=\"$I.png\" border=\"1\" width=\"400\">" >> $html
   echo "$I<br></td>" >> $html

   N=$(($N+1))

   if [ $((N%$ncol)) -eq 0 ]; then
      echo "</tr><tr>" >> $html
   fi
done

echo -e "</tr></table>" >> $html
echo -e "</body>\n</html>" >> $html

echo "... finished."

finishing the story of dratted nitrobenzene

A few days ago I finally finished writing the manuscript for the paper on nitrobenzene photochemistry. I will post a link here as soon as it is published and I dont want to spoil to much. However, a few of the most remarkable outcomes I want to share:

1. If you are planning on doing quantum chemical calculations on molecules including a nitro-group, don't use CC2 for obtaining ground state geometries, it will give you very poor results. I think in general, one should not use CC2 for the ground-state just because it is a coupled cluster approach. Obviously, MP2 is cheaper and seemingly also better/more stable. Even some of the developers of CC2 think so! In my experience, B2PLYP is even better, in particular in complicated electronic situations at pretty much the same or less cost than MP2 (The KS-SCF usually converges much faster and more stable than HF-SCF, which is true in particular for triplet states). For all singly excited-states of nitrobenzene, the excitation energies of TDA-B2PLYP are throughout very close to those obtained with ADC3. Eventually, I am really impressed by the performance of this double-hybrid functional, which it delivers at essentially no cost compared to the expensive ADC3.
2. After doing A LOT of calculations with every theoretical methods I could get my hands on, I think I finally understood what the problem of nitrobenzene is: According to MOM-CCSD(T)/cc-pVTZ, the second lowest excited singlet state of nitrobenzene (1 B1, npi>pi*) has large double excitation character of 50 % according to the ADC3 state composition. Since this state is closely related to the lowest triplet state for which the problems persist, its still unclear which is the lowest triplet state of nitrobenzene. All I can say after doing all these calculations (ADC3, MOM-CCSD(T), EOM-CCSD, DFT-MRCI, CAS-NEVPT2(14/11), TDA-B2PLYP) is that the lowest two triplet states (one n>pi* and one npi>pi*) are very close together regarding their vertical excitation energies (n>pi ~ 0.2 eV above npi>pi) and if geometrical relaxation is considered, they become virtually degenerate (dE < 0.05 eV).
3. I found an experimental spectrum of nitrobenzene vapor published in 1964 [S. Nagakura, M. Kojima, Y. Maruyama, J. Mol. Spec. 13 p. 174], which shows a peak exactly at the excitation energy of 4.4 eV predicted by MOM-CCSD(T) for the doubly excited npi>pi* singlet state. On the basis of a comparison to benzene, Nagakura et Al assigned this peak to an other state. According to my calculations (ADC3 and EOM-CCSD), however, this other state (2 B1) neither has any oscillator strength (~ 0.0002) nor does the energy fit (4.7 eV instead of 4.4 eV). Furthermore is the vibrational splitting of this peak in bare benzene 930 /cm just in the same ballpark compared to the 860 /cm they found for this peak in the spectrum of nitrobenzene. So I calculated the vibrational splitting for the doubly excited 1 B1 state and voi-la: 864 /cm (B2PLYP), which fits much better. After all, I'm pretty sure the first visible peak in the spectrum of nitrobenzene vapor is due to a doubly excited state, which is really unexpected, I guess.

 So far, so good. Now I can turn back to implementing C-PCM for the ADC module of Q-Chem and I already know, on which molecule I will try it out :)

continuing the story of dratted nitrobenzene

ω(r[NO])

I finished investigating the influence of the NO bond length (symmetric stretch) onto excitation energies at ADC(3)/SVP level of theory! The result is as expected: due to the population of the π* orbital at the nitro group (see first post concerning nitrobenzene) in all lower excited-states, the influence of this bond length is quite large in the following way: along the symmetric NO stretching mode, the energy of the ground state largely increases, while the absolute energies of the excited states decrease (see Figure 1).

Figure 1) Run of the absolute ground- and excited state energies along symmetric NO-stretching coordinate starting from the CCSD/def2-TZVP optimized structure. The scan was a rigid one, i.e. without reoptimizing the structural paramaters for each r(NO) but keeping them at the value of the original geometry.

For the vertical excitation energies (ω) these effects add up and the ω of all states drop rapidly by about 0.1 eV/pm (Figure 2), which explains for the large difference between CC2 geometry with 125 pm NO-bonds and the CCSD geometry with 120 pm NO-bonds.

Figure 2) Run of the vertical excitation energies along the symmetric NO-stretching coordinate. The gradient is about 0.1 eV/pm. The equilibrium value for the NO distance is 120 pm at the CCSD/def2-TZVP level of theory (compare Figure 1). One can clearly distinguish between the local excited state of the nitro group and the valence excited states, since the latter have a slightly smaller gradient.

 Now I know why and how strong excitation energies depend on the NO bond length, but one very important question remains: Why is the ADC(3)  result for the excitation energy  of the lowest excited singlet state (almost 4.1 eV) so far off from the experiment (3.65 eV)? From a third order method like ADC(3) one would expect 0.1-0.3 eV deviation at max, in particular for singly excited states. Maybe the reason is lying in the dynamics of the molecule, cause in the experiment, this bond vibrates like any other bond. I will have a look into this using born-oppenheimer ab-initio molecular dynamic in the next post.

 so long!

Accidential Accuracy

 Why chocolate makes you fat but TNT doesn't

TNB: an explosive that is
closely related to TNT

Yesterday we carried out some (as we thought) very approximate calculations for a little 2 hrs quantum chemistry workshop that we will give during the annual meeting of my graduate school next week.

I thought it would be fun to calculate and compare detonation energy of an explosive to the energy content of "chocolate", which basically is sugar and fat 50:50. In doing so, one will find that you could never get fat from TNT because it actually contains only a small fraction of the chemical energy that is stored in chocolate.

a mono sugar: Glucose
(sucrose is too big)

Since we only have two hours for the whole workshop we've chosen to employ an approximate level of theory from the density functional theory (DFT) corner called B3LYP with a small SVP basis on model compounds. Furthermore, we decided to ignore all entropic, thermal and environmental influences to keep it simple. For the detonation of TNB we used a simplified reaction equation yielding only CO, N2 and H2 which circumvents the calculation of electronically complex (radical) nitroxides and such. At this level of theory, all calculations take about 1 hour on a modern laptop computer.

a fat model: Octanoic acid
(real triglycerides are to big)

The trickiest part to get right is the energy of liquid water, which is needed since we have physiological conditions. Here, we extrapolated from a few calculations (one water, one water using a solvent model (PCM), two water + PCM, four water+PCM, eight water+PCM) to to energy of liquid water.

Despite the crude model and level of theory we hoped to achieve at least qualitative agreement with the experiment, but then something not completely unexpected happened.

For comparison (taken from wikipedia)
- The experimental detonation energy of TNB is 3.8-3.9 kJ/g
- The physiological energy content of sugar is 17 kJ/g
- that of fat is 39 kJ/g.

Now look at the B3LYP/SVP results:

results of the B3LYP/SVP thermochemistry calculations

I think the explanation for this astounding accuracy is again a fortuitous cancellation of errors: The intermolecular interactions we neglect on the left side of the equation are obviously of the same or very similar size as the entropic and thermal corrections we neglect on the right side of the equation.

After all, its just spooky how well error cancellation can work out, in particular for DFT/B3LYP. This is not always a good thing! In my last post on nitrobenzene for example, you can read how error-cancellation can be quite problematic.

dratted nitrobenzene!

isosurfaces of the difference densities (DD) and molecular orbitals (MOs) of the lowest electronic states of nitrobenzene

this simple-looking but nasty system had me thinking a lot for the past months.

the problem with its theoretical investigation in a single sentence goes as follows: even high-level and usually very reliable methods (e.g. perturbative coupled-cluster of second order, CC2) systematically deliver largely inaccurate results for very simply properties such as the ground state geometry and vertical excitation energies.

Therefore, expensive third order methods such as the Algebraic Diagrammatic Construction scheme of third order ADC(3) and (Equation of Motion) Coupled-Cluster Singles and Doubles (EOM)-CCSD are required for the vertical excitation energies and the ground- and excited state geometries, respectively.

It took me almost a year (I wasn't constantly working on the topic) to realize that there's something wrong with CC2 due to something that is quite typical for the field I'm working in: an elaborate and fortuitous cancellation of errors. In case of NB and CC2 it goes like that:
If you calculate the molecular ground state geometry of NB with CC2 and subsequently use this geometry to calculate vertical excitation energies with linear-response (LR-)CC2, the outcome is just fine. Fits experimental values without any large unexpected errors and everything seems fine and consistent.
However, if you compare the CC2 geometry to experimental measurements or results from CCSD calculations, you will find that one essential parameter it is systematically wrong (N-O bonds are 125 instead of 122 pm) and it turns out: The vertical excitation energies are only in agreement with the experiment for exactly this wrong geometry.
If experimental parameters or a CCSD optimized geometry is used, the excitation energies with CC2 are all 0.3-0.5 eV too high. Now that I know this,  I started to systematically investigate the influence of geometry, which threw up more questions than it answered.

In the next posts I will probably discuss the influence of geometry on the vertical excitation energies and explain all the shiny pictures above.

So long!