Category Archives: Models

Natural Transition Orbitals (NTOs) Gaussian


The canonical molecular orbital depiction of an electronic transition is often a messy business in terms of a ‘chemical‘ interpretation of ‘which electrons‘ go from ‘which occupied orbitals‘ to ‘which virtual orbitals‘.

Natural Transition Orbitals provide a more intuitive picture of the orbitals, whether mixed or not, involved in any hole-particle excitation. This transformation is particularly useful when working with the excited states of molecules with extensively delocalized chromophores or multiple chromophoric sites. The elegance of the NTO method relies on its simplicity: separate unitary transformations are performed on the occupied and on the virtual set of orbitals in order to get a localized picture of the transition density matrix.

[1] R. L. Martin, J. Chem. Phys., 2003, DOI:10.1063/1.1558471.

In Gaussian09:
After running a TD-DFT calculation with the keyword TD(Nstates=n) (where n = number of states to be requested) we need to take that result and launch a new calculation for the NTOs but lets take it one step at a time. As an example here’s phenylalanine which was already optimized to a minimum at the B3LYP/6-31G(d,p) level of theory. If we take that geometry and launch a new calculation with the TD(Nstates=40) in the route section we obtain the UV-Vis spectra and the output looks like this (only the first three states are shown):

Excitation energies and oscillator strengths:

Excited State 1: Singlet-A 5.3875 eV 230.13 nm f=0.0015 <S**2>=0.000
42 -> 46 0.17123
42 -> 47 0.12277
43 -> 46 -0.40383
44 -> 45 0.50838
44 -> 47 0.11008
This state for optimization and/or second-order correction.
Total Energy, E(TD-HF/TD-KS) = -554.614073682
Copying the excited state density for this state as the 1-particle RhoCI density.

Excited State 2: Singlet-A 5.5137 eV 224.86 nm f=0.0138 <S**2>=0.000
41 -> 45 -0.20800
41 -> 47 0.24015
42 -> 45 0.32656
42 -> 46 0.10906
42 -> 47 -0.24401
43 -> 45 0.20598
43 -> 47 -0.14839
44 -> 45 -0.15344
44 -> 47 0.34182

Excited State 3: Singlet-A 5.9254 eV 209.24 nm f=0.0042 <S**2>=0.000
41 -> 45 0.11844
41 -> 47 -0.12539
42 -> 45 -0.10401
42 -> 47 0.16068
43 -> 45 -0.27532
43 -> 46 -0.11640
43 -> 47 0.16780
44 -> 45 -0.18555
44 -> 46 -0.29184
44 -> 47 0.43124

The oscillator strength is listed on each Excited State as “f” and it is a measure of the probability of that excitation to occur. If we look at the third one for this phenylalanine we see f=0.0042, a very low probability, but aside from that the following list shows what orbital transitions compose that excitation and with what energy, so the first line indicates a transition from orbital 41 (HOMO-3) to orbital 45 (LUMO); there are 10 such transitions composing that excitation, visualizing them all with canonical orbitals is not an intuitive picture, so lets try the NTO approach, we’re going to take excitation #10 for phenylalanine as an example just because it has a higher oscillation strength:

%chk=Excited State 10: Singlet-A 7.1048 eV 174.51 nm f=0.3651 <S**2>=0.000
41 -> 45 0.35347
41 -> 47 0.34685
42 -> 45 0.10215
42 -> 46 0.17248
42 -> 47 0.13523
43 -> 45 -0.26596
43 -> 47 -0.22995
44 -> 46 0.23277

Each set of NTOs for each transition must be calculated separately. First, copy you filename.chk file from the TD-DFT result to a new one and name it after the Nth state of interest as shown below (state 10 in this case). NOTE: In the route section, replace N with the number of the excitation of interest according to the results in filename.log. Run separately for each transition your interested in:

#chk=state10.chk

#p B3LYP/6-31G(d,p) Geom=AllCheck Guess=(Read,Only) Density=(Check,Transition=N) Pop=(Minimal,NTO,SaveNTO)

0 1
--blank line--

By requesting SaveNTO, the canonical orbitals in the state10.chk file are replaced with the NTOs for the 10th excitation, this makes it easier to plot since most visualizers just plot whatever set of orbitals they read in the chk file but if they find the canonical MOs then one would need to do some re-processing of them. This is much more straightforward.

Now we format our chk files into fchk with the formchk utility:

formchk -3 filename.chk filename.fchk
formchk -3 state10.chk state10.fchk

If we open filename.fchk (the file where the original TD-DFT calculation is located) with GaussView we can plot all orbitals involved in excited state number ten, those would be seven orbitals from 41 (HOMO-3) to 47 (LUMO+2) as shown in figure 1.

Figure 1. Canonical orbitals involved in the 10th excited state according to the TD-DFT calculation

If we now open state10.fchk we see that the numbers at the side of the orbitals are not their energy but their occupation number particular to this state of interest, so we only need to plot those with highest occupations, in our example those are orbitals 44 and 45 (HOMO and LUMO) which have occupations = 0.81186; you may include 43 and 46 (HOMO-1 and LUMO+1, respectively) for a much more complete description (occupations = 0.18223) but we’re still dealing with 4 orbitals instead of 7.

Figure 2. Natural Transition Orbitals for Phenylalanine. Orbital 44 (particle) and Orbital 45 (hole) exhibit the largest occupations for Excited State No. 10

The NTO transition 44 -> 45 is far easier to conceptualize than all the 10 combinations given in the canonical basis from the direct TD-DFT calculation. TD-DFT provides us with the correct transitions, NTOs just paint us a picture more readily available to the chemist mindset.

NOTE: for G09 revC and above, the %OldChk option is available, I haven’t personally tried it but using it to specify where the excitations are located and then write the NTOs of interest into a new chk file in the following way, thus eliminating the need of copying the original chk file for each state:

%OldChk=filename.chk
%chk=stateN.chk

NTOs are based on the Natural Hybrid orbitals vision by Löwdin and others, and it is said to be so straightforward that it has been re-discovered from time to time. Be that as it may, the NTO visualization provides a much clearer vision of the excitations occurring during a TD calculation.

Thanks for reading, stay home and stay safe during these harsh days everyone. Please share, rate and comment this and other posts.

DFT Estimation of pKb Values – New Paper in JCIM


As a continuation of our previous work on estimating pKa values from DFT calculations for carboxylic acids, we now present the complementary pKb values for amino groups by the same method, and the coupling of both methodologies for predicting the isoelectric point -pI- values of amino acids as a proof of concept.

Analogously to our work on pKa, we now used the Minimum Surface Electrostatic Potentia, VS,min, as a descriptor of the availability of Nitrogen’s lone pair and correlated it with the experimental basicity of a large number of amines, separated into three groups: primary, secondary and tertiary amines.

Interestingly, the correlation coefficient between experimental and calculated pKb values decreases in the following order: primary (R2 = 0.9519) > secondary (R2 = 0.9112) > tertiary (R2 = 0.8172). This could be due to steric effects, the change in s-character of the lone pair or just plain old selection bias. Nevertheless, there is a good correlation between both values and the resulting equations can predict the pKb value of an amino group within less of a unit, which is very good for a statistical method that does not require the calculation of a full thermodynamic cycle.

We then took thirteen amino acids (those without titratable side chains) and calculated simultaneously VS,min and VS,max for the amino and the carboxyl group (this latter with the use of equation 2 from our previous work published in Molecules MDPI) and the arithmetical average of both gave us their corresponding pI values with an agreement of less than one unit.

This work is now available at the Journal of Chemical Information and Modeling (DOI: 10.1021/acs.jcim.9b01173); as always a shoutout is due to the people working on it: Leonardo “Leo” Lugo, Gustavo “Gus” Mondragón and leading the charge Dr. Jacinto Sandoval-Lira.

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)

eq1 (1)

If you’re trying to calculate interaction energies, remember to perform these corrections for every component in your supramolecular assembly (eq. 2)

eq2(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.

Figure1

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.

 

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:

2017-08-09 20.29.45

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:

ζ1a12s + b12px + c12py + d12pz

ζ2a22s + b22px + c22py + d22pz

ζ3a32s + b32px + c32py + d32pz

ζ4a42s + 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/√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.

No, seriously, why can’t orbitals be observed?


The concept of electronic orbital has become such a useful and engraved tool in understanding chemical structure and reactivity that it has almost become one of those things whose original meaning has been lost and replaced for a utilitarian concept, one which is not bad in itself but that may lead to some wrong conclusions when certain fundamental facts are overlooked.

Last week a wrote -what I thought was- a humorous post on this topic because a couple of weeks ago a viewpoint in JPC-A was published by Pham and Gordon on the possibility of observing molecular orbitals through microscopy methods, which elicited a ‘seriously? again?‘ reaction from me, since I distinctly remember the Nature article by Zuo from the year 2000 when I just had entered graduate school. The article is titled “direct observation of d-orbital holes.” We discussed this paper in class and the discussion it prompted was very interesting at various levels: for starters, the allegedly observed d-orbital was strikingly similar to a dz2, which we had learned in class (thanks, prof. Carlos Amador!) that is actually a linear combination of d(z2-x2) and d(z2-y2) orbitals, a mathematical -lets say- trick to conform to spectroscopic observations.

Pham and Gordon are pretty clear in their first paragraph: “The wave function amplitude Ψ*Ψ is interpreted as the probability density. All observable atomic or molecular properties are determined by the probability and a corresponding quantum mechanical operator, not by the wave function itself. Wave functions, even exact wave functions, are not observables.” There is even another problem, about which I wrote a post long time ago: orbitals are non-unique, this means that I could get a set of orbitals by solving the Schrödinger equation for any given molecule and then perform a unit transformation on them (such as renormalizing them, re-orthonormalizing them to get a localized version, or even hybridizing them) and the electronic density derived from them would be the same! In quantum mechanical terms this means that the probability density associated with the wave function internal product, Ψ*Ψ, is not changed upon unit transformations; why then would a specific version be “observed” under a microscope? As Pham and Gordon state more eloquently it has to do with the Density of States (DOS) rather than with the orbitals. Furthermore, an orbital, or more precisely a spinorbital, is conveniently (in math terms) separated into a radial, an angular and a spin component R(r)Ylm(θ,φ)σ(α,β) with the angular part given by the spherical harmonic functions Ylm(θ,φ), which in turn -when plotted in spherical coordinates- create the famous lobes we all chemists know and love. Zuo’s observation claim was based on the resemblance of the observed density to the angular part of an atomic orbital. Another thing, orbitals have phases, no experimental observation claims to have resolved those.

Now, I may be entering a dangerous comparison but, can you observe a 2? If you say you just did, well, that “2” is just a symbol used to represent a quantity: two, the cardinality of a set containing two elements. You might as well depict such quantity as “II” or “⋅⋅” but still cannot observe “a two”. (If any mathematician is reading this, please, be gentle.) I know a number and a function are different, sorry if I’m just rambling here and overextending a metaphor.

Pretending to having observed an orbital through direct experimental methods is to neglect the Born interpretation of the wave function, Heisenberg’s uncertainty principle and even Schrödinger’s cat! (I know, I know, Schrödinger came up with this gedankenexperiment in order to refute the Copenhagen interpretation of quantum mechanics, but it seems like after all the cat is still not out of the box!)

So, the take home message from the viewpoint in JPC is that molecular properties are defined by the expected values of a given wave function for a specific quantum mechanical operator of the property under investigation and not from the wave function itself. Wave functions are not observables and although some imaging techniques seem to accomplish a formidable task the physical impossibility hints to a misinterpretation of facts.

I think I’ll write more about this in a future post but for now, my take home message is to keep in mind that orbitals are wave functions and therefore are not more observable (as in imaging) than a partition function is in statistical mechanics.

The “art” of finding Transition States Part 2


Last week we posted some insights on finding Transitions States in Gaussian 09 in order to evaluate a given reaction mechanism. A stepwise methodology is tried to achieve and this time we’ll wrap the post with two flow charts trying to synthesize the information given. It must be stressed that knowledge about the chemistry of the reaction is of paramount importance since G09 cannot guess the structure connecting two minima on its own but rather needs our help from our chemical intuition. So, without further ado here is the remainder of Guillermo’s post.

 

METHOD 3. QST3. For this method, you provide the coordinates of your reagents, products and TS (in that order) and G09 uses the QST3 method to find the first order saddle point. As for QST2 the numbering scheme must match for all the atoms in your three sets of coordinates, again, use the connection editor to verify it. Here is an example of the input file.

link 0
--blank line--
#p b3lyp/6-31G(d,p) opt=(qst3,calcfc) geom=connectivity freq=noraman
--blank line--
Charge Multiplicity
Coordinates of reagents
--blank line—
Charge Multiplicity
Coordinates of products
--blank line--
Charge Multiplicity
Coordinates of TS
--blank line---

As I previously mentioned, it happens that you find a first order saddle point but does not correspond to the TS you want, you find an imaginary vibration that is not the one for the bond you are forming or breaking. For these cases, I suggest you to take that TS structure and manually modify the region that is causing you trouble, then use method 2.

METHOD 4. When the previous methods fail to yield your desired TS, the brute force way is to acquire the potential energy surface (PES) and visually locate your possible TS. The task is to perform a rigid PES scan, for this, the molecular structure must be defined using z-matrix. Here is an example of the input file.

link 0
--blank line--
#p b3lyp/6-31G(d,p) scan test geom=connectivity
--blank line--
Charge Multiplicity
Z-matrix of reagents (or products)
--blank line--

In the Z-matrix section you must specify which variables (B, A or D) you want to modify. First, locate the variables you want to modify (distance B, angle A, or dihedral angle D). Then modify those lines within the Z-matrix, here is an example.

B1       1.41     3          0.05
A1       104.5   2          1.0

What you are specifying with this is that the variable B1 (a distance) is going to be stepped 3 times by 0.05. Then variable A1 (an angle) is going to be stepped 2 times by 1.0. Thus, a total of 12 energy evaluations will be performed. At the end of the calculation open the .log file in gaussview and in Results choose the Scan… option. This will open a 3D surface where you should locate the saddle point, this is an educated guess, so take the structure you think corresponds to your TS and use it for method 2.

I have not fully explored this method so I encourage you to go to Gaussian.com and thoroughly review it.

Once you have found your TS structure and via the imaginary vibration confirmed that is the one you are looking for the next step is to verify that your TS connects both your reagents and products in the potential energy surface. For this, an Intrinsic Reaction Coordinate (IRC) calculation must be performed. Here is an example of the input file for the IRC.

link 0
--blank line--
#p b3lyp/6-31G(d,p) irc=calcfc geom=connectivity
--blank line--
Charge Multiplicity
Coordinates of TS
--blank line--

With this input, you ask for an IRC calculation, the default numbers of steps are 20 for each side of your TS in the PES; you must specify the coordinates of your TS or take them from the .chk file of your optimization. In addition, an initial force constant calculation must be made. It often occurs that the calculation fails in the correction step, thus, for complicated cases I hardly suggest to use irc=calcall, this will consume very long time (even days) but there is a 95% guaranty. If the number of points is insufficient you can put more within the route section, here is such an example for a complicated case.

link 0
--blank line--
#p b3lyp/6-31G(d,p) irc=(calcall,maxpoints=80) geom=connectivity
--blank line--
Charge Multiplicity
Coordinates of TS
--blank line--

With this route section, you are asking to perform an IRC calculation with 80 points on each side of the PES, calculating the force constants at every point. For an even complicated case try adding the scf=qc keyword in the route section, quadratic convergence often works better for IRC calculations.

The ‘art’ of finding Transition States Part 1


Guillermo CaballeroGuillermo Caballero, a graduate student from this lab, has written this two-part post on the nuances to be considered when searching for transition states in the theoretical assessment of reaction mechanisms. He’s been quite successful in getting beautiful energy profiles for organic reaction mechanisms, some of which have even explained why some reactions do not occur! A paper in Tetrahedron has just been accepted but we’ll talk about it in another post. I wanted Guillermo to share his insight into this hard practice of computational chemistry so he wrote the following post. Enjoy!

 

Yes, finding a transition state (TS) can be one of the most challenging tasks in computational chemistry, it requires both a good choice of keywords in your route section and all of your chemical intuition as well. Herein I give you some good tricks when you have to find a transition state using Gaussian 09 Rev. D1

METHOD 1. The first option you should try is to use the opt=qst2 keyword. With this method you provide the structures of your reagents and your products, then the program uses the quadratic synchronous transit algorithm to find a possible transition state structure and then optimize it to a first order saddle point. Here is an example of the input file.

link 0
--blank line--
#p b3lyp/6-31G(d,p) opt=qst2 geom=connectivity freq=noraman
--blank line--
Charge Multiplicity
Coordinates of reagents
--blank line--
Charge Multiplicity
Coordinates of products
--blank line---

It is mandatory that the numbering must be the same in the reagents and the products otherwise the calculation will crash. To verify that the label for a given atom is the same in reagents and products you can go to Edit, then Connection. This opens a new window were you can manually modify the numbering scheme. I suggest you to work in a split window in gaussview so you can see at the same time your reagents and products.

The keyword freq=noraman is used to calculate the frequencies for your optimized structure, it is important because for a TS you must only observe one imaginary frequency, if not, then that is not a TS and you have to use another method. It also occurs that despite you find a first order saddle point, the imaginary frequency does not correspond to the bond forming or bond breaking in your TS, thus, you should use another method. I will give you advice later in the text for when this happens. When you use the noraman in this keyword you are not calculating the Raman frequencies, which for the purpose of a TS is unnecessary and saves computing time. Frequency analysis MUST be performed AT THE VERY SAME LEVEL OF THEORY at which the optimization is performed.

The main advantage for using the qst2 option is that if your calculation is going to crash, it generally crashes at the beginning, in the moment of guessing your transition state structure. Once the program have a guess, it starts the optimization. I suggest you to ask the algorithm to calculate the force constants once, this generally improves on the convergence, it will take slightly more time depending on the size of your structure but it pays off. The keyword in the route section is opt=(qst2,calcfc). Indeed, I hardly encourage you to use the calcfc keyword in any optimization you want to run.

METHOD 2. If method 1 does not work, my next advice is to use the opt=ts keyword. For this method, the coordinates in your input file are those for the TS structure. Here is an example of the input file.

link 0
--blank line--
#p b3lyp/6-31G(d,p) opt=ts geom=connectivity freq=noraman
--blank line--
Charge Multiplicity
Coordinates of TS
--blank line--

The question that arises here is how should I get the coordinates for my TS? Well, honestly this is not a trivial task, here is where you use all the chemistry you know. For example, you can start with the coordinates of your reagents and manually get them closer. If you are forming a bond whose length is to be 1.5Å, then I suggest you to have that length in 1.6Å in your TS. Sometimes this becomes trial and error but the most accurate your TS structure is, based on your chemical knowledge, the easiest to find your TS will be. As another example, if you want to find a TS for a [1,5]-sigmatropic reaction a good TS structure will be putting the hydrogen atom that migrates in the middle point through the way. I have to insist, this method hardly depends on your imagination to elucidate a TS and on your chemistry background.

Most of the time when you use the opt=ts keyword the calculations crashes because of an error in the number of eigenvalues, you can avoid it adding noeigen to the route section; here is an example of the input file, I encourage you to use this method.

link 0
--blank line--
#p b3lyp/6-31G(d,p) opt=(ts,noeigen,calcfc) geom=connectivity freq=noraman
--blank line--
Charge Multiplicity
Coordinates of TS
--blank line--

If you have problems in the optimization steps I suggest you to ask the algorithm to calculate the force constants in every step of the optimization opt=(ts,noeigen,calcall) this is quite a harsh method because will consume long computing time but works well for small molecules and for complicated TSs to find.

Another ‘tricky’ way to get your coordinates for your TS is to run the qst2 calculation, then if it fails, take the second- or the third-step coordinates and used them as a ‘pre-optimized’ set of coordinates for this method.

By the way, here is another useful trick. If you are evaluating a group of TSs, let’s say, if you are varying a functional group among the group, focus on finding the TS for the simplest case, then use this optimized TS as a template where you add the moieties and use this this method. This works pretty well.

For this post we’ll leave it up to here and post the rest of Guillermo’s tricks and advice on finding TS structures next week when we’ll also discuss the use of IRC calculations and some considerations on energy corrections when plotting the full energy profile. In the mean time please take the time to rate, like and share this and other posts.

Thanks for reading!

 

Quantifying σ-Holes – G09 and MultiWFN


Some atomic properties such as an atomic charge are isotropic, but every now and then some derivations of them become anisotropic, for example the plotting of the Molecular Electrostatic Potential (MEP) on the electron density surface can exhibit some anisotropic behavior; quantifying it can be a bit challenging.

It is well known that halogen atoms such as Chlorine can form so-called halogen-bonds of the type R-Cl-R in crystals with a near perfect 180° angle. This finding has lead to the discovery of σ-holes in halogens. σ-Holes are electrophilic portions of the anisotropic electrostatic potential in an otherwise nucleophilic atom. Recently, Guillermo “Memo” Caballero and I calculated the MEP for a series of trichloromethyl-containing compounds at the MP2/cc-pVQZ level of theory and the mapping shows evidence of such σ-holes as seen in Figure1. Those small blue portions on an otherwise red atom indicate that some electron density is missing in that position, which by the way is located at 180° away from the carbon atom.

But having the picture is not enough. We want to quantify just how strong are those σ-holes to effectively attract a nucleophile and perhaps perform some chemistry on the C-Cl bond. That’s when we resorted to MultiWFN, a Multifunctional Wavefunction Analyzer developed by Tian Lu (卢天) at the Beijing Kein Research Center for Natural Sciences. You can check the project leader list of publications here. Among many other capabilities, MultiWFN is able to print details about properties along a surface.

In order to work with MultiWFN you need to generate a *.wfn file, if you have a previous Gaussian calculation for which you want to analyze their surface you can run a guess=only calculation in order to extract the wavefunction from the checkpoint file. Here is a dummy of the input for such calculation

%chk=oldfile.chk
# output=wfn geom=check guess=(only,read) density=current
--blank--
Title Card
--blank--
0 1
--blank--
filename.wfn
--blank--

In our case, having a post-Hartree-Fock calculation, the use of density=current is mandatory to get the MP2 density matrix and not just the HF one. Running this calculation will generate the file filename.wfn which is now used with MultiWFN. When starting MultiWFN you get to see a terminal window like the one below in which you are asked to input the path of your wfn file:

1

After loading it you will get the following window with the various options available. Type 12// (these two slashes are mandatory) to get the quantitative analysis of molecular surface option.

2

Then you will be asked to define some elements of that surface (we used the default options 0)

3

The following screenshot shows the results section in which several maxima  and minima of electrostatic potential were found (7 and 11 in our case); a star is placed on the side of the global maximum. The value of the MEP at those points is given in Hartrees, eV and kcal/mol which I personally hate because there isn’t such a thing as a mole of ‘potentials’ (same argument as giving an orbital’s energy in kcal/mol, moles of what? orbitals? Personally, I don’t like it even if its valid).4

Their visualizer is activated through the option 0 and although it is far from pretty it is quite good enough to find the numbers corresponding to maxima and minima of the MEP on the isodensity surface. If we look for the maxima then we find for our example (CHCl3) that a maximum is located in front of each Cl atom in a straight line from the C atom. Now we get to put a number on the mapped isosurface provided by Gaussian or even import the file into Chimera.

5ccl3cn

We are still working our way around MultiWFN, I hope we can find the batch option, it would be most useful. In the mean time, Guillermo and I will continue to search for σ-holes in chlorinated reagents. Thanks to Guillermo for his ongoing work in this and other topics within the realm of organic reactivity.

Have you any suggestions or ideas to work with MultiWFN? Please share them!

A personal artistic impression of the CompChem Landscape 


In a nutshell, computational chemistry models are about depicting, reproducing and predicting the electronic-based molecular reality. I had this conversation with my students last week and at some point I drew a parallel between them and art in terms of how such reality is approached.

Semi empirical methods
Prehistoric wall paintings depict a coarse aspect of reality without any detail but nevertheless we can draw some conclusions from the images. In the most sophisticated of these images, the cave paintings in Altamira, we can discern a bison, or could it be a bull? but definitely not a giraffe nor a whale, most in the same way Hückel´s method provides an ad hoc picture of π electron density without any regard of the σ portion of the electron density or the conformational possibilities (s-cis and s-trans 1,3-butadiene have the same Hückel description).

More sophisticated semi-empirical Hamiltonians like PM3 or PM6 have better parametrizations and hence yield better results. We are still replacing a lot of information for experimental or adjusted parameters but we still cannot truly adopt it as truthful. Take this pre-medieval painting of one of the first Kings of England, Aelred the Unready. It is, by today standards, a good children´s drawing and not a royal portrait, we now see more detail and can discern many more features yielding a better description of a human figure than those found in Altamira or Egypt.

Ethelred_the_Unready

Æthelred the Unready King of England ca. 1000 BC

 

Hartree-Fock
HF is the simplest of ab initio methods, meaning that no experimental results or adjustable parameters are introduced. Even more so, from the HF equations for a multi-electron system that complies with Pauli’s exclusion principle the exchange operator arises as a new quantum feature of matter with no classical analogue. Still, there are some shortcomings. Correlation energy is disregarded and most results vary according to the basis set employed. Take the impressionist movement, specially in France: In Monet´s Lady with Umbrella we have a more complicated composition, we observe many more features and although we have a better description of color composition some details, like her face, remain obscure. The impressionists are characterized by their broad strokes, the thicker the strokes the harder it is to observe details similar to what happens in HF when we change from a small to a large basis set, respectively.

Claude_Monet_023

Woman with a parasol – Her face or the identity of the flowers at her feet are indistinguishable yet we might be safe to say its springtime.

CI (Configurations Interaction)
Extension of HF to a multi-reference method yields better results. In CI we take the original guess wavefunction -as expressed through a Slater Determinant- and extend it with one or many more wavefunctions; thus a linear combination of Slater Determinants gives rise to a broader description of the ground state because other electronic configurations are involved to include more details like the ionic and covalent pictures (configurations). The more terms we include the more real the results feel. If we take classical figurative paintings we have a similar result; most of these paintings are constituted of many elements and the more realistically each element is captured the more real the whole composition looks even if some are just merely indicated.

Flaming_June,_by_Frederic_Lord_Leighton_(1830-1896)

Flaming June by Lord Leighton – Extreme details on the fabrics and the sea in the background makes us oblivious to the less detailed foot

CCSD(T) full-CI, CASPT2

In Edwards Much’s the scream, we might think we have lost some information again and went back to impressionism but we know this is actually an expressionist painting; we can now not only observe details of the figurative portion of the image but Munch has captured his subject´s fear in the form of distorsions on the subjective reality. In this way, CCSD(T), full-CI and CASPT2 methods provide a description of the ground as well as the excited states which -in experimental reality- are only accessed through a perturbation of the elecron density by electromagnetic radiation. Something resembling radiation has perturbated the subject in The Scream rendering him frightened and wondering how to return to his ground state or if such thing will be even possible.

The_Scream

The Scream by Edward Munch – what sort of perturbation got this guy’s fears out?

 

Density Functional Methods

At least due to its widespread use, DFT has risen as the preferred method. One of the reasons behind its success is the reduced computing time when compared to previous ab initio methods. So DFT is pretty much like photography, in which reality is captured in full but only apparently after selecting a given lens, an exposition, a filter, shutter speed and the occasional Photoshop for correcting issues such as aliasing. In photography, as in DFT, all details concerning the procedure or method for capturing an uncanny reproduction of reality must be stated in every case for reproduction purposes.

Now, in the end it all comes down to Magritte’s Pipe. Ceci n’est pas une pipe -or, ‘this is not a pipe’- reminds us that painting as with modeling we don’t get reality but rather a depiction of it. In this famous painting we look at an image that in our heads resembles that of a pipe but we cannot grab it, fill it with tobacco and smoke it.

MagrittePipe

This isn’t even Magritte’s painting! Let alone a pipe

The image above is a digital file, which translated becomes a scaled reproduction of an image painted by Magritte in which we see the 2D projection of the image of an object that reminds us of a pipe. In fact, the real name of this work is The Treachery of Images, definitely quite an epistemology problem on perception and knowledge but before I get too metaphysical I should finish this post.

Can you find where cubism or surrealism should be placed? with MPn methods, perhaps?

Comparing the Relative Stability of Non-Equivalent Molecules


How do you compare the stability of two or more compounds which differ in some central atom(s)?

If you simply calculate the energy of both compounds you get a misleading answer since the number of electrons is different from one to the next -in fact, the answer is not so much misleading as it is erroneous. Take compounds 1 and 2 shown in figure 1, for example. Compound 1 was recently synthesized characterized through X-Ray crystallography by my friend Dr. Monica Moya’s group; compound 2 doesn’t exist and we want to know why – or at least know if it is relatively unstable respect to 1.

Figure 1. Compound 1 exists but compound 2 is apparently less stable. Is it?

Figure 1. Compound 1 exists but compound 2 is apparently less stable. Is it?

Although stoichiometry is the same, varying only by the substitution of Ga by Al the number of electrons is quite different. We then made the following assumption: Since the atomic radii of Ga and Al are quite similar (according to the CCDC their respective covalent radii are 122[4] and 121[3] pm), relative stability must rely on the bonding properties rendering 2 harder to obtain, at least through the method used for 1. The total energy for compound 1 was calculated at the M06-2X/6-31G(d,p) level of theory; then both Al atoms were changed by Ga and the total energy was calculated again at the same level. Separately, the energy of isolated Ga and Al atoms were calculated. Compensating the number of electrons was now a simple algebraic problem:

ΔE = E(MnBxOy) – nEM + nEM’ – E(M’nBxOy)

 The absolute energy difference E1 – E2 is staggering due to the excess of 36 electrons in 2. But after this compensation procedure we now have a more reliable result of ΔE value of ca. 81 kcal/mol in favor of compound 1. In strict sense, we performed geometry optimizations at various stages: first on compound 1 to remove the distorsions due to the crystal field and then on the substituted compound 2 to make sure Ga atoms would find a right fit in the molecule but since their covalent radii are similar, no significant changes in the overall geometry were observed confirming the previous assumption.

We now have the value of the energy difference between 1 and 2 and other similar cases, the next step is to find the distal causes of the relative stability which may rely on the bonding properties of the Ga-O bond respect to the Al-O bonds.

What do you think? Is there another method you can share for tackling this problem? Please share your thoughts on the comments section.

Thanks for reading.

%d bloggers like this: