Blog Archives
Calculating NMR shifts – Short and Long Ways
Nuclear Magnetic Resonance is a most powerful tool for elucidating the structure of diamagnetic compounds, which makes it practically universal for the study of organic chemistry, therefore the calculation of 1H and 13C chemical shifts, as well as coupling constants, is extremely helpful in the assignment of measured signals on a spectrum to an actual functional group.
Several packages offer an additive (group contribution) empirical approach to the calculation of chemical shifts (ChemDraw, Isis, ChemSketch, etc.) but they are usually only partially accurate for the simplest molecules and no insight is provided for the more interesting effects of long distance interactions (vide infra) so quantum mechanical calculations are really the way to go.
With Gaussian the calculation is fairly simple just use the NMR keyword in the route section in order to calculate the NMR shielding tensors for relevant nuclei. Bear in mind that an optimized structure with a large basis set is required in order to get the best results, also the use of an implicit solvation model goes a long way. The output displays the value of the total isotropic magnetic shielding for each nucleus in ppm (image taken from the Gaussian website):
Magnetic shielding (ppm): 1 C Isotropic = 57.7345 Anisotropy = 194.4092 XX= 48.4143 YX= .0000 ZX= .0000 XY= .0000 YY= -62.5514 ZY= .0000 XZ= .0000 YZ= .0000 ZZ= 187.3406 2 H Isotropic = 23.9397 Anisotropy = 5.2745 XX= 27.3287 YX= .0000 ZX= .0000 XY= .0000 YY= 24.0670 ZY= .0000 XZ= .0000 YZ= .0000 ZZ= 20.4233
Now, here is why this is the long way; in order for these values to be meaningful they need to be contrasted with a reference, which experimentally for 1H and 13C is tetramethylsilane, TMS. This means you have to perform the same calculation for TMS at -preferably- the same level of theory used for the sample and substract the corresponding values for either H or C accordingly. Only then the chemical shifts will read as something we can all remember from basic analytical chemistry class.
GaussView 6.0 provides a shortcut; open the Results menu, select NMR and in the new window there is a dropdown menu for selecting the nucleus and a second menu for selecting a reference. In the case of hydrogen the available references are TMS calculated with the HF and B3LYP methods. The SCF – GIAO plot will show the assignments to each atom, the integration simulation and a reference curve if desired.
The chemical shifts obtained this far will be a good approximation and will allow you to assign any peaks in any given spectrum but still not be completely accurate though. The reasons behind the numerical deviations from calculated and experimental values are many, from the chosen method to solvent interactions or basis set limitations, scaling factors are needed; that’s when you can ask the Cheshire Cat which way to go
If you don’t know where you are going any road will get you there.
Lewis Carroll – Alice in Wonderland
Well, not really. The Chemical Shift Repository for computed NMR scaling factors, with Coupling Constants Added Too (aka CHESHIRE CCAT) provides with straight directions on how to correct your computed NMR chemical shifts according to the level of theory without the need to calculate the NMR shielding tensor for the reference compound (usually TMS as pointed out earlier). In a nutshell, the group of Prof. Dean Tantillo (UC Davis) has collected a large number of isotropic magnetic shielding values and plotted them against experimental chemical shifts. Just go to their scaling factors page and check all their linear regressions and use the values that more closely approach to your needs, there are also all kinds of scripts and spreadsheets to make your job even easier. Of course, if you make use of their website don’t forget to give the proper credit by including these references in your paper.
We’ve recently published an interesting study in which the 1H – 19F coupling constants were calculated via the long way (I was just recently made aware of CHESHIRE CCAT by Dr. Jacinto Sandoval who knows all kinds of web resources for computational chemistry calculations) as well as their conformational dependence for some substituted 2-aza-carbazoles (fig. 1).
Journal of Molecular Structure Vol 1176, 15 January 2019, Pages 562-566
The paper is published in the Journal of Molecular Structure. In this study we used the GIAO NMR computations to assign the peaks on an otherwise cluttered spectrum in which the signals were overlapping due to conformational variations arising from the rotation of the C-C bond which re-orients the F atoms in the fluorophenyl grou from the H atom in the carbazole. After the calculations and the scans were made assigning the peaks became a straightforward task even without the use of scaling factors. We are now expanding these calculations to more complex systems and will contrast both methods in this space. Stay tuned.
Post Calculation Addition of Empirical Dispersion – Fixing interaction energies
Calculation of interaction energies is one of those things people are more concerned with and is also something mostly done wrong. The so called ‘gold standard‘ according to Pavel Hobza for calculating supramolecular interaction energies is the CCSD(T)/CBS level of theory, which is highly impractical for most cases beyond 50 or so light atoms. Basis set extrapolation methods and inclusion of electronic correlation with MP2 methods yield excellent results but they are not nonetheless almost as time consuming as CC. DFT methods in general are terrible and still are the most widely used tools for electronic structure calculations due to their competitive computing times and the wide availability of schemes for including terms which help describe various kinds of interactions. The most important ingredients needed to get a decent to good interaction energies values calculated with DFT methods are correlation and dispersion. The first part can be recreated by a good correlation functional and the use of empirical dispersion takes care of the latter shortcoming, dramatically improving the results for interaction energies even for lousy functionals such as the infamous B3LYP. The results still wont be of benchmark quality but still the deviations from the gold standard will be shortened significantly, thus becoming more quantitatively reliable.
There is an online tool for calculating and adding the empirical dispersion from Grimme’s group to a calculation which originally lacked it. In the link below you can upload your calculation, select the basis set and functionals employed originally in it, the desired damping model and you get in return the corrected energy through a geometrical-Counterpoise correction and Grimme’s empirical dispersion function, D3, of which I have previously written here.
The gCP-D3 Webservice is located at: http://wwwtc.thch.uni-bonn.de/
The platform is entirely straightforward to use and it works with xyz, turbomole, orca and gaussian output files. The concept is very simple, a both gCP and D3 contributions are computed in the selected basis set and added to the uncorrected DFT (or HF) energy (eq. 1)
(1)
If you’re trying to calculate interaction energies, remember to perform these corrections for every component in your supramolecular assembly (eq. 2)
(2)
Here’s a screen capture of the outcome after uploading a G09 log file for the simplest of options B3LYP/6-31G(d), a decomposed energy is shown at the left while a 3D interactive Jmol rendering of your molecule is shown at the right. Also, various links to the literature explaining the details of these calculations are available in the top menu.
I’m currently writing a book chapter on methods for calculating ineraction energies so expect many more posts like this. A special mention to Dr. Jacinto Sandoval, who is working with us as a postdoc researcher, for bringing this platform to my attention, I was apparently living under a rock.
DFT Textbook in Spanish by Dr. José Cerón-Carrasco
Today’s science is published mostly in English, which means that non-English speakers must first tackle the language barrier before sharing their scientific ideas and results with the community; this blog is a proof that non-native-English speakers such as myself cannot outreach a large audience in another language.
For young scientists learning English is a must nowadays but it shouldn’t shy students away from learning science in their own native tongues. To that end, the noble effort by Dr. José Cerón-Carrasco from Universidad Católica San Antonio de Murcia, in Spain, of writing a DFT textbook in Spanish constitutes a remarkable resource for Spanish-speaking computational chemistry students because it is not only a clear and concise introduction to ab initio and DFT methods but because it was also self published and written directly in Spanish. His book “Introducción a los métodos DFT: Descifrando B3LYP sin morir en el intento” is now available in Amazon. Dr. Cerón-Carrasco was very kind to invite me to write a prologue for his book, I’m very thankful to him for this opportunity.
Así que para los estudiantes hispanoparlantes hay ahora un muy valioso recurso para aprender DFT sin morir en el intento gracias al esfuerzo y la mente del Dr. José Pedro Cerón Carrasco a quien le agradezco haberme compartido la primicia de su libro
¡Salud y olé!
Another Great Year at the Lab! 2017
2017 was a complicated year for various reasons here in Mexico (and some personal health issues) but nonetheless I’m very proud of the performance of everyone at the lab whose hard work and great skills keep pushing our research forward.
Four new members joined the team and have presented their work at the national meeting for CompChem for the first time. Also, for the first time, one of my students, Gustavo Mondragón, gave a talk at this meeting with great success about his research on the Fenna Matthews Olson complex of photosynthetic bacteria.
The opportunity to attend WATOC at Munich presented me the great chance to meet wonderful people from around the world and was even kindly and undeservingly invited to write the prologue for an introductory DFT book by Prof. Pedro Cerón from Spain. I hope to Jeep up with the collaborations abroad such as the one with the Mirkin group at Nortgwestern and the one with my dear friend Kunsagi-Mate Sándor at Pecsi Tudomanyegyetem (Hungary), among many others; I’m thankful for their trust in our capabilities.
Two members got their BSc degrees, Marco an Durbis, the latter also single handedly paved the way for us to develop a new research line on the in silico drug developing front; his relentless work has also been praised by the QSAR team at the Institute of Chemistry with which he has collaborated by performing toxicity calculations for the agrochemical industry as well as by designing educational courses aimed to the dissemination of our work and QSAR in general among regulatory offices and potential clients. We’re sad to see him go next fall but at the same time we’re glad to know his scientific skills will further develop.
I cannot thank the team enough: Alejandra Barrera, Gustavo Mondragón, Durbis Castillo, Fernando Uribe, Juan Guzman, Alberto Olmedo, Eduardo Cruz, Ricardo Loaiza and Marco Garcia; may 2018 be a great year for all of you.
And to all the readers thank you for your kind words, I’m glad this little space which is about to become nine years old is regarded as useful; to all of you I wish a great 2018!
-
- three UBPs
Python scripts for calculating Fukui Indexes
One of the most popular posts in this blog has to do with calculating Fukui indexes, however, when dealing with a large number of molecules, our described methodology can become cumbersome since it requires to manually extract the population analysis from two or three different output files and then performing the arithmetic on them separately with a spreadsheet or something.
Our new team member Ricardo Loaiza has written a python script that takes the three aforementioned files and yields a .csv file with the calculated Fukui indexes, and it even points out which of the atoms exhibit the largest values so if you have a large molecule you don’t have to manually check for them. We have also a batch version which takes all the files in any given directory and performs the Fukui calculations for each, provided it can find file triads with the naming requirements described below.
Output files must be named filename.log (the N electrons reference state), filename_plus.log (the state with N+1 electrons) and filename_minus.log (the N-1 electrons state). Another restriction is that so far these scripts only work with NBO population analysis as provided by the NBO3.1 program available in the various versions of Gaussian. I imagine the listing is similar in NBO5.x and NBO6.x and so it should work if you do the population analysis with them.
The syntax for the single molecule version is:
python fukui.py filename.log filename_minus.log filename_plus.log
For the batch version is:
./fukuiPorLote.sh
(Por Lote means In Batch in Spanish.)
These scripts are available via GitHub. We hope you find them useful, and you do please let us know whether here at the comments section or at our GitHub site.
XVI Mexican Meeting on Phys.Chem.
A yearly tradition of this Comp.Chem. lab and many others throughout our nation is to attend the Mexican Meeting on Theoretical Physical Chemistry to share news, progress and also a few drinks and laughs. This year the RMFQT was held in Puebla and although unfortunately I was not able to attend this lab was proudly represented by its current members. Gustavo Mondragón gave a talk about his progress on his photosynthesis research linking to the previous work of María Eugenia Sandoval already presented in previous editions; kudos to Gustavo for performing remarkably and thanks to all those who gave us their valuable feedback and criticism. Also, five posters were presented successfully, I can only thank the entire team for representing our laboratory in such an admirable way, and a special mention to the junior members, I hope this was the first of many scientific events they attend and may you deeply enjoy each one of them.
Among the invited speakers, the RMFQT had the honor to welcome Prof. John Perdew (yes, the P in PBE); the team took the opportunity of getting a lovely picture with him.
Here is the official presentation of the newest members of our group:
Alejandra Barrera (hyperpolarizabilty calculations on hypothetical poly-calyx[n]arenes for the search of NLO materials)
Fernando Uribe (Interaction energy calculations for non-canonical nucleotides)
Juan Guzmán (Reaction mechanisms calculations for catalyzed organic reactions)
We thank the organizing committee for giving us the opportunity to actively participate in this edition of the RMFQT, we eagerly await for next year as every year.
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!
All you wanted to know about Hybrid Orbitals…
… but were afraid to ask
or
How I learned to stop worrying and not caring that much about hybridization.
The math behind orbital hybridization is fairly simple as I’ll try to show below, but first let me give my praise once again to the formidable Linus Pauling, whose creation of this model built a bridge between quantum mechanics and chemistry; I often say Pauling was the first Quantum Chemist (Gilbert N. Lewis’ fans, please settle down). Hybrid orbitals are therefore a way to create a basis that better suits the geometry formed by the bonds around a given atom and not the result of a process in which atomic orbitals transform themselves for better sterical fitting, or like I’ve said before, the C atom in CH4 is sp3 hybridized because CH4 is tetrahedral and not the other way around. Jack Simmons put it better in his book:
Taken from “Quantum Mechanics in Chemistry” by Jack Simmons
The atomic orbitals we all know and love are the set of solutions to the Schrödinger equation for the Hydrogen atom and more generally they are solutions to the hydrogen-like atoms for which the value of Z in the potential term of the Hamiltonian changes according to each element’s atomic number.
Since the Hamiltonian, and any other quantum mechanical operator for that matter, is a Hermitian operator, any given linear combination of wave functions that are solutions to it, will also be an acceptable solution. Therefore, since the 2s and 2p valence orbitals of Carbon do not point towards the edges of a tetrahedron they don’t offer a suitable basis for explaining the geometry of methane; even more so these atomic orbitals are not degenerate and there is no reason to assume all C-H bonds in methane aren’t equal. However we can come up with a linear combination of them that might and at the same time will be a solution to the Schrödinger equation of the hydrogen-like atom.
Ok, so we need four degenerate orbitals which we’ll name ζi and formulate them as linear combinations of the C atom valence orbitals:
ζ1= a12s + b12px + c12py + d12pz
ζ2= a22s + b22px + c22py + d22pz
ζ3= a32s + b32px + c32py + d32pz
ζ4= a42s + b42px + c42py + d42pz
to comply with equivalency lets set a1 = a2 = a3 = a4 and normalize them:
a12 + a22 + a32 + a42 = 1 ∴ ai = 1/√4
Lets take ζ1 to be directed along the z axis so b1 = c1 = 0
ζ1 = 1/√4(2s) + d12pz
since ζ1 must be normalized the sum of the squares of the coefficients is equal to 1:
1/4 + d12 = 1;
d1 = √3/2
Therefore the first hybrid orbital looks like:
ζ1 = 1/√4(2s) +√3/2(2pz)
We now set the second hybrid orbital on the xz plane, therefore c2 = 0
ζ2 = 1/√4(2s) + b22px + d22pz
since these hybrid orbitals must comply with all the conditions of atomic orbitals they should also be orthonormal:
〈ζ1|ζ2〉 = δ1,2 = 0
1/4 + d2√3/2 = 0
d2 = –1/2√3
our second hybrid orbital is almost complete, we are only missing the value of b2:
ζ2 = 1/√4(2s) +b22px +-1/2√3(2pz)
again we make use of the normalization condition:
1/4 + b22 + 1/12 = 1; b2 = √2/√3
Finally, our second hybrid orbital takes the following form:
ζ2 = 1/√4(2s) +√2/√3(2px) –1/√12(2pz)
The procedure to obtain the remaining two hybrid orbitals is the same but I’d like to stop here and analyze the relative direction ζ1 and ζ2 take from each other. To that end, we take the angular part of the hydrogen-like atomic orbitals involved in the linear combinations we just found. Let us remember the canonical form of atomic orbitals and explicitly show the spherical harmonic functions to which the 2s, 2px, and 2pz atomic orbitals correspond:
ψ2s = (1/4π)½R(r)
ψ2px = (3/4π)½sinθcosφR(r)
ψ2pz = (3/4π)½cosθR(r)
we substitute these in ζ2 and factorize R(r) and 1/√(4π)
ζ2 = (R(r)/√(4π))[1/√4 + √2 sinθcosφ –√3/√12cosθ]
We differentiate ζ2 respect to θ, and set it to zero to find the maximum value of θ respect to the z axis we get the angle between the first to hybrid orbitals ζ1 and ζ2 (remember that ζ1 is projected entirely over the z axis)
dζ2/dθ = (R(r)/√(4π))[√2 cosθ –√3/√12sinθ] = 0
sinθ/cosθ = tanθ = -√8
θ = -70.53°,
but since θ is measured from the z axis towards the xy plane this result is equivalent to the complementary angle 180.0° – 70.53° = 109.47° which is exactly the angle between the C-H bonds in methane we all know! and we didn’t need to invoke the unpairing of electrons in full orbitals, their promotion of any electron into empty orbitals nor the ‘reorganization‘ of said orbitals into new ones. Orbital hybridization is nothing but a mathematical tool to find a set of orbitals which comply with the experimental observation and that is the important thing here!
To summarize, you can take any number of orbitals and build any linear combination you want, in order to comply with the observed geometry. Furthermore, no matter what hybridization scheme you follow, you still take the entire orbital, you cannot take half of it because they are basis functions. That is why you should never believe that any atom exhibits something like an sp2.5 hybridization just because their bond angles lie between 109 and 120°. Take a vector v = xi+yj+zk, even if you specify it to be v = 1/2i that means x = 1/2, not that you took half of the unit vector i, and it doesn’t mean you took nothing of j and k but rather than y = z = 0.
This was a very lengthy post so please let me know if you read it all the way through by commenting, liking, or sharing. Thanks for reading.
