Mittwoch, 22. Januar 2014

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:

Next, the user reads the geometry (e.g. from any xyz file) into vmd (vmd, 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:
# 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


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


if [ $1 -eq 2 ]
  isov2=`echo $isov'/8' | bc -l`
#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 ]
  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
   echo "Please enter 2 or 3 for # of surfaces!"
   echo "Falling back to 3 surfaces with standard values"
  echo "Using standard isovalues: "
  echo "Values: "$isov $isov2 $isov3
  echo "material change opacity Glass3 0.400000" > $out

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

for I in *$ifmt
   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


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

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

echo "... finished."

Montag, 20. Januar 2014

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 :)