excited-state solvent-field equilibration in pratice

Sketch of a solvent-field equilibration for
the cationic ethene dimer.

In my last post, I introduced the cationic ethene dimer shown on the right to illustrate and evaluate a solvent model for excited states. Despite its neat size and simplicity,  quantum-chemical calculations on this system were a bit tricky to converge onto the desired states. This is due to its open-shell (doublet) electronic structure exhibiting multiple, energetically almost degenerate SCF solutions. While these problems can be circumvented for the gas-phase calculation with one little trick, the actual solvent-field equilibration requires for more drastic measures, some of which I want to present in this post. Here are the three main problems and how I approached them:

1. The SCF converges on an undesired solution: Without a PCM as well as for low dielectric constants and large basis set, the initial SCF converges onto the symmetrically charged (both ethenes + ½) solution, instead of the asymmetrically charged one I want to investigate.
   
   To convince the SCF algorithm to converge onto the asymmetrically charged solution, I use sequential jobs. In my favorite quantum-chemistry program, this can be triggered by adding a line “@@@” after the first input and appending the second one. The advantage of this is that results stored on disk, like e.g. MOs, state-densities and PCM surface-charges are available in the following step. I want to use the first job to create a set of asymmetric orbitals, which I employ as initial guess in the second one. For this purpose, I shorten the C-H bonds of one of the molecules by 0.1 A to introduce an asymmetry that triggers the SCF to converge onto the asymmetrical solution. Using the resulting orbitals as starting point for the SCF, I can give the second job the undistorted geometry (otherwise it would by default use the one from the first job). Starting from the asymmetrical guess orbitals, the second SCF also converges onto the desired solution, even for the now symmetric geometry. This trick was necessary for all gas-phase calculations, as well as for the calculations with large (cc-pVQZ) or augmented (aug-cc-pVDZ) basis sets, even if a PCM with high epsilon is used.
   
   
2. The energetic ordering of the excited states switches: While the charge-inverse excited state is one of the higher lying states in the ground-state equilibrated solvent field (2.-4., depending on epsilon and the basis), it becomes the lowest one as soon as its solvent field is relaxed.
   
   Changes in the energetic ordering of the excited-states during excited state geometry optimizations are very common, if not inherent to the problem and I guess every theoretical photochemist knows the case. Using symmetry in the calculation can help, but does not necessarily solve the problem. Surprisingly few quantum-chemistry codes feature a decent root homing algorithm to identify excited-states via their overlap with the results from the previous step to guide the optimization. For the solvent-field optimizations of this predictable toy system the solution is rather trivial: Since I have not yet implemented a proper iteration loop, but use a script to generate sequential input files with a fixed number of iteration steps, I can just adapt the state_to_opt variable separately in each of the iteration steps. Nevertheless, I wanted to point out how practical root homing would be as a common feature in quantum-chemistry codes.
   
   
3. The SCF "follows" the solvent field: The first SCF calculation in the solvent-field of the excited state converges onto a solution resembling the former excited-state wavefunction, since the latter is much lower in energy in the new, excited-state relaxed solvent field.
   
   This one was the trickiest one to tackle. Although I do per default use the orbitals of the previous step as a guess for the solvent-field iterations all of the available SCF algorithms slowly crawl to the solution of the energetically much lower state. Obviously, the excited-state solvent field is just too attractive for the positive charge. In my favorite quantum-chemistry program, however, there is a very nice feature which I want to introduce in a little more detail in the following. The maximum overlap method (MOM) is basically a kind of inter-step root-homing for the SCF, which was developed by Andrew Gilbert and Peter Gill for the calculation of excited-states wavefunctions. For this purpose, converged HF orbitals from a previous calculation are rotated in such a way that they resemble the excited state, e.g. via  "excitation" of an electron from the HOMO to the LUMO. One drawback of the method is that for this to succeed, the overlap between the ground- and respective excited state wavefunction has to be rather small. Hence, the approach does in general only work for dark excited states. If you are interested any further: I discussed the advantages and drawbacks of this method for the description of excited states in combination with coupled-cluster theory somewhat more extensively in a work on nitrobenzene. For the problem at hand, i.e. forcing the SCF to converge onto the local minimum provided in the guess orbitals, MOM works like a charm and I consider making its use the default for my equilibrium solvation approach.

After all, the results I obtained with these tricks are promising: If the solvent-fields are fully relaxed for the respective, correlated MP2 or ADC2 density, the energy difference between the ground- and charge-inverse excited state stay within 0.01 eV of the gas-phase value for the whole range of epsilon:

Energy difference between the solvent equilibrated ground (MP2/cc-pVTZ) and charge-inverse excited (ADC2/cc-pVTZ) states (left y-axis) as a function of the dielectric constant (x-axis). For clarity, I have subtracted the gas-phase value of 0.34 eV from the results. For relation, the total solute-solvent interaction energy for one of the states is given with respect to the right y-axis.

Since this post has already become quite long and the day quite old I will continue with a more detailed survey of the results as well as a description of the differences between the PTE, PTD and PTED schemes in my next post.

So long!
Jan

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.