Blog Archives

The Evolution of Photosynthesis

Recently, the journal ACS Central Science asked me to write a viewpoint for their First Reactions section about a research article by Prof. Alán Aspuru-Guzik from Harvard University on the evolution of the Fenna-Matthews-Olson (FMO) complex. It was a very rewarding experience to write this piece since we are very close to having our own work on FMO published as well (stay tuned!). The FMO complex remains a great research opportunity for understanding photosynthesis and thus the origin of life itself.

In said article, Aspuru-Guzik’s team climbed their way up a computationally generated phylogenetic tree for the FMO from different green sulfur bacteria by creating small successive mutations on the protein at a time while also calculating their photochemical properties. The idea is pretty simple and brilliant: perform a series of “educated guesses” on the structure of FMO’s ancestors (there are no fossil records of FMO so this ‘educated guesses’ are the next best thing) and find at what point the photochemistry goes awry. In the end the question is which led the way? did the photochemistry led the way of the evolution of FMO or did the evolution of FMO led to improved photochemistry?

Since both the article and viewpoint are both published as open access by the ACS, I wont take too much space here re-writing the whole thing and will instead exhort you to read them both.

Thanks for doing so!


Chemistry Makes the Chemical

The compound shown below in figure 1 is listed by Aldrich as 4,5,6,7-tetrahydroindole, but is it really?


Fig 1. An indole?

To a hardcore organic chemist it is clear that this is not an indole but a pyrrole because  the lack of aromaticity in the fused ring gives this molecule the same reactivity as 2,3-diethyl pyrrole.  If you search the ChemSpider database for ‘tetrahydroindole’ the search returns the following compound with the identical chemical formula C8H11N but with a different hydrogenation pattern: 2,3,3a,4-Tetrahydro-1H-indole


Fig 2. Also listed as an indole

The real indole, upon an electrophilic attack, behaves as a free enamine yielding the product shown in figure 3 in which the substitution occurs in position 3. This compound cannot undergo an Aromatic Electrophilic Susbstitution since that would imply the formation of a sigma complex which would disrupt the aromaticity.


On the contrary, the corresponding pyrrole is substituted in position 2


These differences in reactivity towards electrophiles are easily rationalized when we plot their HOMO orbitals (calculated at the M062X/def2TZVP level of theory):

If we calculate the Fukui indexes at the same level of theory we get the highest value for susceptibility towards an electrophilic attack as follows: 0.20 for C(3) in indole and 0.25 for C(2) in pyrrole, consistent with the previous reaction schemes.

So, why is it listed as an indole? why would anyone search for it under that name? Nobody thinks about cyclohexane as 1,3,5-trihydrobenzene. According to my good friend and colleague Dr. Moisés Romero most names for heterocyles are kept even after such dramatic chemical changes due to historical and mnemonic reasons even when the reactivity is entirely different. This is only a nomenclature issue that we have inherited from the times of Hantzsch more than a century ago. We’ve become used to keeping the trivial (or should I say arbitrary) names and further use them as derivations but this could pose an epistemological problem if students cannot recognize which heterocycle presents which reactivity.

So, in a nutshell:

Chemistry makes the chemical and not the structure.

A thing we all know but sometimes is overlooked for the sake of simplicity.

WATOC 2017

Last week the WATOC congress in Munich was a lot of fun. Our poster on photosynthesis had a great turnout and got a lot of positive feedback as well as many thought provoking questions. One of the highlights of my time there was seeing my former students and knowing they’re all leading successful and happy grad-student lives in Europe, I’m so very proud of them. It was great to connect with old friends and making new ones; a big thank you to all the readers of this little blog who took the time to come and say hi, I’m very glad the blog has been helpful to you.

Better recounts of WATOC 2017 can be found in the great Rzepa’s blog here and here.

Below there is an image of our poster (some typos persist).


See you all in 2020!

Photosynthesis and Singlet Fission – #WATOC2017 PO1-296

If you work in the field of photovoltaics or polyacene photochemistry, then you are probably aware of the Singlet Fission (SF) phenomenon. SF can be broadly described as the process where an excited singlet state decays to a couple of degenerate coupled triplet states (via a multiexcitonic state) with roughly half the energy of the original singlet state, which in principle could be centered in two neighboring molecules; this generates two holes with a single photon, i.e. twice the current albeit at half the voltage (Fig 1).


Jablonski’s Diagram for SF

It could also be viewed as the inverse process to triplet-triplet annihilation. An important requirement for SF is that the two triplets to which the singlet decays must be coupled in a 1(TT) state, otherwise the process is spin-forbidden. Unfortunately (from a computational perspective) this also means that the 3(TT) and 5(TT) states are present and should be taken into account, and when it comes to chlorophyll derivatives the task quickly scales.

SF has been observed in polyacenes but so far the only photosynthetic pigments that have proven to exhibit SF are some carotene derivatives; so what about chlorophyll derivatives? For a -very- long time now, we have explored the possibility of finding a naturally-occurring, chlorophyll-based, photosynthetic system in which SF could be possible.

But first things first; The methodology: It was soon enough clear, from María Eugenia Sandoval’s MSc thesis, that TD-DFT wasn’t going to be enough to capture the whole description of the coupled states which give rise to SF. It was then that we started our collaboration with SF expert, Prof. David Casanova from the Basque Country University at Donostia, who suggested the use of Restricted Active Space – Spin Flip in order to account properly for the spin change during decay of the singlet excited state. A set of optimized bacteriochlorophyll-a molecules (BChl-a) were oriented ad-hoc so their Qy transition dipole moments were either parallel or perpendicular; the rate to which SF could be in principle present yielded that both molecules should be in a parallel Qy dipole moments configuration. When translated to a naturally-occurring system we sought in two systems: The Fenna-Matthews-Olson complex (FMO) containing 7 BChl-a molecules and a chlorosome from a mutant photosynthetic bacteria made up of 600 Bchl-d molecules (Fig 2). The FMO complex is a trimeric pigment-protein complex which lies between the antennae complex and the reaction center in green sulfur dependent photosynthetic bacteria such as P. aestuarii or C. tepidium, serving thus as a molecular wire in which is known that the excitonic transfer occurs with quantum coherence, i.e. virtually no energy loss which led us to believe SF could be an operating mechanism. So far it seems it is not present. However, for a crystallographic BChl-d dimer present in the chlorosome it could actually occur even when in competition with fluorescence.


FMO Complex. Trimer (left), monomer (center), pigments (right)


BChQRU chlorosome. 600 Bchl-d molecules

I will keep on blogging more -numerical and computational- details about these results and hopefully about its publication but for now I will wrap this post by giving credit where credit is due: This whole project has been tackled by our former lab member María Eugenia “Maru” Sandoval and Gustavo Mondragón. Finally, after much struggle, we are presenting our results at WATOC 2017 next week on Monday 28th at poster session 01 (PO1-296), so please stop by to say hi and comment on our work so we can improve it and bring it home!

Some Comp.Chem. Tweeps

Out of some +1000 twitter accounts I follow about a quarter are related computational chemistry. The following public list isn’t comprehensive and prone to errors and contains researchers, programmers, students, journals, products and companies who gravitate around the use of in silico methods for the understanding and design of chemical and biochemical compounds.


Dealing with Spin Contamination

Most organic chemistry deals with closed shell calculations, but every once in a while you want to calculate carbenes, free radicals or radical transition states coming from a homolytic bond break, which means your structure is now open shell.

Closed shell systems are characterized by having doubly occupied molecular orbitals, that is to say the calculation is ‘restricted’: Two electrons with opposite spin occupy the same orbital. In open shell systems, unrestricted calculations have a complete set of orbitals for the electrons with alpha spin and another set for those with beta spin. Spin contamination arises from the fact that wavefunctions obtained from unrestricted calculations are no longer eigenfunctions of the total spin operator <S^2>. In other words, one obtains an artificial mixture of spin states; up until now we’re dealing only with single reference methods. With each step of the SCF procedure the value of <S^2> is calculated and compared to s(s+1) where s is half the number of unpaired electrons (0.75 for a radical and 2.0 for triplets, and so on); if a large deviation between these two numbers is found, the then calculation stops.

Gaussian includes an annihilation step during SCF to reduce the amount of spin contamination but it’s not 100% reliable. Spin contaminated wavefunctions aren’t reliable and lead to errors in geometries, energies and population analyses.

One solution to overcome spin contamination is using Restricted Open Shell calculations (ROHF, ROMP2, etc.) for which singly occupied orbitals is used for the unpaired electrons and doubly occupied ones for the rest. These calculations are far more expensive than the unrestricted ones and energies for the unpaired electrons (the interesting ones) are unreliable, specially spin polarization is lost since dynamical correlation is hardly accounted for. The IOP(5/14=2) in Gaussian uses the annihilated wavefunction for the population analysis if acceptable but since Mulliken’s method is not reliable either I don’t advice it anyway. 

The case of DFT is different since rho.alpha and rho.beta can be separated (similarly to the case of unrestricted ab initio calculations), but the fact that both densities are built of Kohn-Sham orbitals and not true canonical orbitals, compensates the contamination somehow. That is not to say that it never shows up in DFT calculations but it is usually less severe, of course for the case of hybrid functional the more HF exchange is included the more important spin contamination may become. 

So, in short, for spin contaminated wavefunctions you want to change from restricted to unrestricted and if that doesn’t work then move to Restricted Open Shell; if using DFT you can use the same scheme and also try changing from hybrid to pure orbitals at the cost of CPU time. There is a last option which is using spin projection methods but I’ll discuss that in a following post. 

Stability of Unnatural DNA – @PCCP #CompChem

As is the case of proteins, the functioning of DNA is highly dependent on its 3D structure and not just only on its sequence but the difference is that protein tertiary structure has an enormous variety whereas DNA is (almost) always a double helix with little variations. The canonical base pairs AT, CG stabilize the famous double helix but the same cannot be guaranteed when non-canonical -unnatural- base pairs (UBPs) are introduced.


Figure 1

When I first took a look at Romesberg’s UBPS, d5SICS and dNaM (throughout the study referred to as X and Y see Fig.1) it was evident that they could not form hydrogen bonds, in the end they’re substituted naphtalenes with no discernible ways of creating a synton like their natural counterparts. That’s when I called Dr. Rodrigo Galindo at Utah University who is one of the developers of the AMBER code and who is very knowledgeable on matters of DNA structure and dynamics; he immediately got on board and soon enough we were launching molecular dynamics simulations and quantum mechanical calculations. That was more than two years ago.

Our latest paper in Phys.Chem.Chem.Phys. deals with the dynamical and structural stability of a DNA strand in which Romesberg’s UBPs are introduced sequentially one pair at a time into Dickerson’s dodecamer (a palindromic sequence) from the Protein Data Bank. Therein d5SICS-dNaM pair were inserted right in the middle forming a trisdecamer; as expected, +10 microseconds molecular dynamics simulations exhibited the same stability as the control dodecamer (Fig.2 left). We didn’t need to go far enough into the substitutions to get the double helix to go awry within a couple of microseconds: Three non-consecutive inclusions of UBPs were enough to get a less regular structure (Fig. 2 right); with five, a globular structure was obtained for which is not possible to get a proper average of the most populated structures.

X and Y don’t form hydrogen bonds so the pairing is pretty much forced by the scaffold of the rest of the DNA’s double helix. There are some controversies as to how X and Y fit together, whether they overlap or just wedge between each other and according to our results, the pairing suggests that a C1-C1′ distance of 11 Å is most stable consistent with the wedging conformation. Still much work is needed to understand the pairing between X and Y and even more so to get a pair of useful UBPs. More papers on this topic in the near future.

I’m putting a new blog out there

As if I didn’t have enough things to do I’m launching a new blog inspired by the #365papers hashtag on Twitter and the blog. In it I’ll hopefully list, write a femto-review of all the papers I read. This new effort is even more daunting than the actual reading of the huge digital pile of papers I have in my Mendeley To-Be-Read folder, the fattest of them all. The papers therein wont be a comprehensive review of Comp.Chem. must-read papers but rather papers relevant to our lab’s research or curiosity.

Maybe I’ll include some papers brought to my attention by the group and they could do the review. The whole endeavor might flop in a few weeks but I want to give it a shot; we’ll see how it mutates and if it survives or not. So far I haven’t managed to review all papers read but maybe this post will prompt to do so if only to save some face. The domain of the new blog is but I think it should have included the word MY at the beginning so as to convey the idea that it is only my own biased reading list. Anyway, if you’re interested share it and subscribe, those post will not be publicized.

Unnatural DNA and Synthetic Biology

Ever since I read the highly praised article by Floyd Romesberg in Nature back in 2013 I got really interested in synthetic biology. In said article, an unnatural base pair (UBP) was not only inserted into a DNA double strand in vivo  but the organism was even able to reproduce the UBPs present in subsequent generations.


Romesberg’s Nucleosides. No Hydrogen bonding is formed between them!

Inserting new unnatural base pairs in DNA works a lot like editing a computer’s code. Inserting a couple UBPs in vitro is like inserting a comment; it wont make a difference but its still there. If the DNA sequence containing the UBPs can be amplified by molecular biology techniques such as PCR it means that a polymerase enzyme is able to recognize it and place it in site, this is equivalent to inserting a ‘hello world’ section into a working code; it will compile but it’s pretty much useless. Inserting these UBPs in vivo means that the organism is able to thrive despite the large deformation in a short section of its genetic code, but having it replicated by the chemical machinery of the nucleus is an amazing feat that only a few molecules could allow.

The ultimate goal of synthetic biology would be to find a UBP which codes effectively and purposefully during translation of DNA.This last feat would be equivalent to inserting a working subroutine in a program with a specific purpose. But not only could the use of UBPs serve for the purposes of expanding the genetic code from a quaternary (base four) to a senary (base six) system: the field of DNA origami could also benefit from having an expansion in the chemical and structural possibilities of the famous double helix; marking and editing a sequence would also become easier by having distinctive sections with nucleotides other than A, T, C and G.

It is precisely in the concept of double helix that our research takes place since the available biochemical machinery for translation and replication can only work on a double helix, else, the repair mechanisms get activated or the DNA will just stop serving its purpose (i.e. the code wont compile).

My good friend, Dr. Rodrigo Galindo and I have worked on the simulation of Romesberg’s UBPs in order to understand the underlying structural, dynamical and electronic causes that made them so successful and to possibly design more efficient UBPs based on a set of general principles. A first paper has been accepted for publication in Phys.Chem.Chem.Phys. and we’re very excited for it; more on that in a future post.

Grimme’s Dispersion DFT-D3 in Gaussian #CompChem

I was just asked if it is possible to perform DFT-D3 calculations in Gaussian and my first answer was to use the following  keyword:


which is available in G16 and G09 only in revision D, apparently.

There are also some overlays that can be used to invoke the use dispersion in various scenarios:

IOp(3/74=x) Exchange and Correlation Potentials





DSD-PBEP86 (double hybrid, DFT-D3).


B2PLYP-D3 (double hybrid, DFT-D3).

B97-D (DFT-D3).

IOp(3/76=x) Mixing of HF and DFT.

-33 PW6B95 and PW6B95-D3 coefficients.

IOp(3/124=x) Empirical dispersion term.




Force dispersion type 3 (Grimme DFT-D3).

Force dispersion type 4 (Grimme DFT-D3(BJ)).

Force dispersion type 5 (Grimme D3, PM7 version).


The D3 correction method of Grimme defines the van der Waals energy like:

$\displaystyle E_{\rm disp} = -\frac{1}{2} \sum_{i=1}^{N_{at}} \sum_{j=1}^{N_{at...
...{6ij}} {r_{ij,{L}}^6} +f_{d,8}(r_{ij,L})\,\frac{C_{8ij}} {r_{ij,L}^8} \right ),$

where coefficients $ C_{6ij}$ are adjusted depending on the geometry of atoms i and j. The damping D3 function for is:

$\displaystyle f_{d,n}(r_{ij}) = \frac{s_n}{1+6(r_{ij}/(s_{R,n}R_{0ij}))^{-\alpha_{n}}},$

where the values of s are adjustable parameters fit for the exchange-correlation functionals used in each calculation.

%d bloggers like this: