excited-state solvent models and a toy system

My latest projects are focused around a solvent model I extended and interfaced to the excited state methods developed in our group during my PhD. But before I write about the toy system I came up with to test my latest work, let me briefly introduce you to polarizable-continuum solvation models. PCMs basically compute a set of point charges on the molecular surface which depend on the molecular electrostatic potential (ESP) to mimic the interaction with solvent molecules. These (point) charges are added to the one-electron Hamiltonian (just like the nuclei) during the self-consistent field (SCF) calculation and updated in every iteration. Eventually, one obtains a wavefunction of a molecule in solution and the interaction between the point charges and the molecular ESP corresponds to the solute-solvent interaction. The single parameter used to describe the solvent (more parameters are used, e.g. for the cavity construction) is its macroscopic dielectric constant epsilon.

Visualization of the PCM surface-charges (blue -, red +) obtained for my beloved example nitrobenzene in its electronic ground state. As one would expect, the surface-charge close to the negatively charged nitro group (front right) is positive - just like the potential nitrobenzene would exhibit in a polar solvent.

 After finishing and extensively evaluating a so-called non-equilibrium (NEq) solvent model for vertical excitation energies against a set of experimental data for nitroaromatics (find it here), I have recently been working on the respective equilibrium variant. In a nutshell, non-equilibrium and equilibrium solvent models can be rationalized in analogy to the Franck-Condon (FC) principle:

• The non-equilibrium case applies to vertical (i.e. very fast) processes like absorption and fluorescence. To be in accordance with the FC principle, you want to clamp the solvents nuclei, while relaxing its electronic degrees of freedom. Within a PCM, this is realized by using the optical dielectric constant (n²) instead of epsilon for a recomputation of the point charges w.r.t. the excited-state ESP. n² is nothing but the macroscopic polarizability ofa solvent at the frequency of visible light, whereas epsilon is the total polarizability (including a rearrangement of the solvent nuclei).
• The equilibrium case applies to long-lived states, i.e. whenever you expect your systems nuclei to relax w.r.t. an excited state potential energy surface. Consequently, you should also relax the solvents nuclei, which is equivalent to recompute the surface charge for the new molecular charge distribution using the full dielectric constant just like for the ground state.

Schematic visualization of equilibrium and non-equilibirum solvation for the example of an excited-state proton transfer in quinolines. The proton transfer can be regarded as the molecules geometrical relaxation, which is accompanied by the solvent equilibration.

Although the explanation and also the equations for the non-equilibrium case are slightly more complicated due to the separation into fast (n²) and slow (epsilon) parts of the polarization, they are much simpler to implement. The reason being that n² of typical solvents is rather small one hence, one can apply perturbation theory an turn it into a correction for excitation energies. For details, see my paper.

The equilibrium case looks straightforward at first glance, but turns out much more complicated w.r.t. its implementation because of the way our excited-state methods work: We compute the Hartree-Fock (HF) ground-state wavefunction (in the form of molecular orbitals, short MOs) and subsequently use it as a basis for the calculation of excited-state wavefunctions and excitation energies with ADC (Algebraic-Diagrammatic Construction for the Polarization Propagator, a configuration-interaction type excited-state method).
Since the HF MOs contain the interaction with PCM surface charges just like they contain the interaction with the nuclei, one needs the PCM surface charges of the excited state before you do the initial HF calculation. Since this is not possible, one has to iterate between HF and ADC until the solvent field and the excited-state ESP are converged, like so:

1. Standard HF calculation with PCM solvent model => MOs
2. ADC calculation using HF MOs => excited-state density
3. HF calculation in the "frozen" surface-charges for ES density=> new MOs
4. ADC calculation with new MOs to obtain new excited-state density
5. repeat from 3 until surface charges/energy are converged 

Eventually, I got the implementation working and the first numbers were looking promising. To check if the code does indeed work correctly, I came up with a simple toy system in which the solvated ground state is identical to an equilibrium solvated charge-transfer excited state: an asymmetrically charged, cationic ethene dimer with a distance of 1 nm:

The cationic ethene dimer as test system. The influence from the PCM is indicated by the boldness of the cavity. As soon as the iterations of the ES reaction field are converged, the system is in the same situation as in the beginning. Also, the non equilibrium situations with charge-tranfer (CT) state the in the field of the ground state (abbreviated  "0" in the scheme) (Nr. 2) and the ground state in the field of the CT state (converged Nr. 3) should be very comparable.

There is only a one practical problem: During the solvent-equilibration iterations (described above), the charge in the ground state jumps back and forth, causing the solvent-equilibration procedure to oscillate instead of converging. Obviously, the HF-SCF calculation in the frozen surface charges (point 3) of the charge-tranfer state converges to the energetically much more favorable solution with the positive charge sitting the ethene in the negatively charged cavity (No surprise, this is exactly what the SCF algorithm is supposed to do). Employing the MOs of the previous iteration step as starting point in the subsequent SCF cycles is not enough do do the trick, the other solution seems too low energetically.

How I was nevertheless able to converge the calculations using a number of tricks will be presented along with some more details on the system and preliminary results in my next post very soon.

So long!
Jan

publishing in JPC

Publishing can be kind of annoying. In particular if journals outsource the last quanta of work required to convert a manuscript prepared with latex into a journal article in PDF format to the authors.
JPC recently decided to require the authors to provide references in manuscript with article titles in title case as well as some other changes. To save at least some authors from the pain of  doing this by hand or searching the internet for the tools, here they are:
1) Install the achemso package and use:
\bibliographystyle{achemso}
\AtBeginDocument{\nocite{achemso-control}}
2) Add these (hopefully self-explanatory) lines in the referenced .bib file:
@Control{achemso-control,
  ctrl-Article-Title  = "yes",
  ctrl-etal-number    = "10",
  ctrl-max-authors    = "10",
}
3) Use the perl script provided and described here:
http://www.stat.berkeley.edu/~paciorek/computingTips/Change_case_your_journal_ti.html
to change all the articles titles in the .bib file to title case.
4) Write an angry mail to ACS and complain.
So long!

Happy New Years!

Apparently, the past few months it was pretty quiet here. The reason was that I was completely occupied writing up a paper on the central project of my Ph.D. project followed by my thesis. Furthermore, there were the formalities of handing it in and organizing the defence, which I definitely underestimated.
After all, I got everything aligned and there will probably more activity in the next time. To begin with, I present to you my take on Theoretical Chemistry, which constitutes the foreword of my thesis.

In the next posts, as soon as I figure out a comfortable way of translating Tex code into HTML, I will publish my take on the Hartree-Fock self-consistent field method, Configuration-Interaction for ground and excited states, Perturbation Theory applied to HF as well as to CI, yielding the Intermediate-State Representation for the description of correlated excited-states.

So long.

On Theoretical Chemistry

 Chemistry is defined as the science of the interaction and interconversion of complex atomic systems called molecules, including metals and salts. Since molecules are the subunits of any common matter, chemistry is ubiquitous and has a long history as scientific discipline. In the early 20th century, it was discovered that atoms themselves are not indivisible as their Greek-descending name (a-tomos, in-divisible) suggests, but are composed of protons, neutrons and electrons. Nowadays, only the latter are still seen as indivisible, or in other words fundamental with respect to the standard model of particle physics. Protons and neutrons, in contrast consist of even smaller subunits: so-called up and down quarks.

These quarks, which  are again fundamental particles, carry fractions of the elementary charge of either plus two thirds (up quarks) or minus one third (down quarks). Protons consist of two up quarks and one down quark and thus carry one positive elementary charge, while for neutrons the charges of one up quark and two down quarks cancel out. The very number of quarks that constitute each proton and neutron does, in combination with their fractional charge, lead to an important coincidence: any atomic nucleus consisting of an arbitrary number of protons and neutrons will carry an exact multiple of the elementary charge and can in turn be neutralized by adding the respective number of electrons. Formally, this process yields the neutral atoms that constitute the periodic table of the elements.

Due to the huge binding energy of quarks, which is the result of the so-called strong interaction between them, they can only be observed in groups. Similarly, with the exception of hydrogen, neutrons and protons are fused together in the atomic core, where they are held in place by the very same short-ranged forces of the strong interaction. Despite the 10^5 times smaller size compared to the atom including the electrons, the nucleus contains virtually the complete mass of an atom. The biggest part of this mass is due to the huge binding energy according to Einstein's equation m = Ec^(-2), whereas the resting mass contributes less than 1%.

Ultimately, atoms may be rationalized as the frozen interaction energy of their elementary building blocks carrying a positive charge, which incidentally is a multiple of the elementary charge. Moreover, this charge exclusively determines which chemical element a nucleus belongs to, and gives rise to the Coulomb potential that attracts the much lighter, negatively charged electrons. It is the delicate interplay of many electrons in the field of complex nuclear arrangements which guides the nuclear motion and in turn the whole of chemistry with its vast number of stable molecules and reactions. Following this line of thought, chemistry could formally be categorized as the sub-field of particle physics concerned with the quantum mechanics of the electromagnetic interaction in complex systems composed of electrons and nuclei.
However, such a categorization would neither reflect the vast number of thinkable and/or stable molecules emerging from the interplay of electrons in the field of the nuclei, nor the great relevance this discipline had long before the connections to the underlying physical equations were discovered. The sheer infinite possibilities to combine atoms to molecules, investigate their properties in creative ways and find empirical solutions to chemical problems (e.g. acid-base models, Lewis-structures gave rise to an independent scientific discipline. Nevertheless, the utilization of the laws of quantum-mechanics to establish a first-principles approach for chemical problems is and has been very promising ever since these connections were first discovered at the beginning of the 20th century, as evident from the following quote:


"The underlying physical laws necessary for the mathematical theory of a large part of physics and the whole of chemistry are thus completely known, and the difficulty is only that the exact application of these laws leads to equations much too complicated to be soluble. It therefore becomes desirable that approximate practical methods of applying quantum mechanics should be developed, which can lead to an explanation of the main features of complex atomic systems without too much computation."

Paul A. M. Dirac, 1929

In the past decades ever faster computers in combination with powerful algorithms and wisely chosen approximations made it possible to solve these insoluble equations on a regular basis. Eventually, the empirical top-down models developed by chemists over hundreds of years could for the first time be challenged by a first-principles bottom-up approach providing a new quality of insights into the elementary steps of chemical reactions such electron transfer as well as into the electronic structure of matter. This first-principles approach is what characterizes the field of theoretical chemistry, which emerged from a fruitful collaboration of chemists, physicists, mathematicians and computer scientists over the past hundred years.

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