Blog Archives

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!

Partial Optimizations with Gaussian09


Sometimes you just need to optimize some fragment or moiety of your molecule for a number of reasons -whether because of its size, your current interest, or to skew the progress of a previous optimization- or maybe you want just some kind of atoms to have their positions optimized. I usually optimize hydrogen atoms when working with crystallographic files but that for some reason I want to preserve the rest of the molecule as refined, in order to keep it under a crystalline field of sorts.
Asking Gaussian to optimize some of the atoms in your molecule requires you to make a list albeit the logic behind it is not quite straightforward to me. This list is invoked by the ReadOptimize keyword in the route section and it includes all atoms by default, you can then further tell G09 which atoms are to be included or excluded from the optimization.
So, for example you want to optimize all atoms EXCEPT hydrogens, then your input should bear the ReadOptimize keyword in the route section and then, at the end of the molecule specification, the following line:

atoms notatoms=H

If you wish to selectively add some atoms to the list while excluding others, here’s an example:

atoms=C H S notatoms=5-8

This list adds, and therefore optimizes, all carbon, hydrogen and sulfur atoms, except atoms 5, 6, 7 and 8, should they be any of the previous elements in the C H S list.
The way I selectively optimize hydrogen atoms is by erasing all atoms from the list -using the noatoms instruction- and then selecting which are to be included in the list -with atoms=H-, but I haven’t tried it with only selecting hydrogen atoms from the start, as in atoms=H

noatoms atoms=H

I probably get very confused because I learned to do this with the now obsolete ReadFreeze keyword; now it sometimes may seem to me like I’m using double negatives or something – please do not optimize all atoms except if they are hydrogen atoms. You can include numbers, ranks or symbols in this list as a final line of your input file.

Common errors (by common I mean I’ve got them):

Lets look at the end of an input I just was working with:

>  AtmSel:  Line=”P  0″
>  Maximum list size exceeded in AddBin.
>  Error termination via Lnk1e in…

AtmSel is the routine which reads the atoms list and I was using a pseudopotential on phosphorous atoms, I placed the atoms list at the end of the file but it should be placed right after the coordinates and the connectivity matrix, should there be one, and thus before any external basis set or pseudopotential or any other specification to be read by Gaussian.

As a sort of test you can use the instruction:

%kjob l103
%chk=myfile.chk
...

at the Link0 section (where your checkpoint is defined). This will kill the job after the link 103 is finished, thus you will only get a list of what parameters were frozen and which were active. Then, if things look ok, you can run the job without the %kjob l103 instruction and get it done.

As usual I hope this helps. Thanks for reading except to those who didn’t read it except for the parts they did read.

Simulation of Raman Spectroscopy and crystal cell effects – Selenium Carboxylate Eur. J. Inorg. Chem.


Computing spectroscopic features of molecules is always an interesting challenge, specially when intermolecular contacts are into play. Take vibrational spectroscopy for instance, all the non-covalent interactions present in a solid will have an important effect on the the calculated frequencies and their intensities. However calculating the spectroscopical properties of a solid quickly becomes a daunting task.

My colleague and friend Dr. Vojtech Jancik asked me to calculate the Raman frequencies for a new compound: Selenoyl bis-carboxylate, which according to him was very hard to obtain due to the very nature of selenium. So we performed various calculations on the isolated molecule to reproduce the measured Raman spectrum but we soon realized that a calculation on the crystal cell was needed if we wanted to get a more thorough picture of the experiment.

The level of theory used was PBEPBE/LANL2DZ. Optimization of the title structure pointed to a low coordination capacity by carboxylate groups as evidenced by the longer Se -O-C=O distances and reduced Wiberg bond indexes. A blue shift was observed for all bands and so we calculated the Raman frequencies at the crystal structure which gave us a better correspondence between spectra. Finally we computed the Raman spectra for the full unit cell comprised of four molecules with which an excellent agreement was obtained (a scaling factor of 0.8 was used).

Unfortunately we failed to further extend this calculation to a larger system with four unit cells and 32 molecules apparently due to insufficient memory; the calculation just stalled and stopped without error after consuming its time in the queue. I’ll try to take a look into it some day.

https://vine.co/v/eaQOH57nADO/embed/simple?audio=1https://platform.vine.co/static/scripts/embed.js

You can read the whole story in: Synthesis and Crystal Structure of the First Selenonyl Bis(carboxylate) SeO2(O2CCH3)2
Lukas Richtera · Vojtech Jancik · Joaquín Barroso‐Flores · Petr Nykel · Jiri Touzin · Jan Taraba

European Journal of Inorganic Chemistry 06/2015; 2015(18):2923–2927. DOI:10.1002/ejic.201500271
https://platform.vine.co/static/scripts/embed.js

Thanks for reading!

Efficient use of the clipboard on GaussView


Editing large molecules on a seemingly simple visualizer as GaussView can be a bit daunting. I’m working on a follow up of that project we recently published in JACS but now we require to attach two macrocycles to the organometallic moiety; the only caveat is that this time we don’t have any crystallographic data with which to start. Generating a 3D model of this structure is already hard enough and even when you managed to do it there are many degrees of freedom that in some cases can lead to unrealistic geometries after optimization.

I recently came across a simple way to edit a large complicated molecule by optimizing the fragments separately and then joining them in a new molecule by using the clipboard. This rather simple method, that I for one had never exploited has just saved me a few good hours.

Copy a molecule (CTRL+C) and it will go to the clipboard as a molecular fragment for which you can define a new hot atom and thus bind it to the other fragment as you would with the regular builder. I strongly suggest to use a “New Molecule Group” instead of editing over an existing molecule. Also, if you are using the “paste” button, observe that it has three different options; you may want to use the last one “append to existing molecule” or you will have your fragments in different windows.

And remember, dihedral angles are your best friends.

Atoms in Molecules (QTAIM) – Flash lesson


As far as population analysis methods goes, the Quantum Theory of Atoms in Molecules (QTAIM) a.k.a Atoms in Molecules (AIM) has become a popular option for defining atomic properties in molecular systems, however, its calculation is a bit tricky and maybe not as straightforward as Mulliken’s or NBO.

Personally I find AIM a philosophical question since, after the introduction of the molecule concept by Stanislao Cannizzaro in 1860 (although previously developed by Amadeo Avogadro who was dead at the time of the Karlsruhe congress), the questions of whether or not an atom retains its identity when bound to others? where does an atom end and the next begins? What are the connections between atoms in a molecule? are truly interesting and far deeper than we usually consider because it takes a big mental leap to think about how matter is organized to give rise to substances. Particularly I’m very interested with the concept of a Molecular Graph which in turn is concerned with the way we “draw lines” to form conceptual molecules. Perhaps in a different post we can go into the detail of the method, which is based in the Laplacian operator of the electron density, but today, I just want to collect the basic steps in getting the most basic AIM answers for any given molecule. Recently, my good friend Pezhman Zarabadi-Poor and I have used rather extensively the following procedure. We hope to have a couple of manuscripts published later on. Therefore, I’ve asked Pezhman to write a sort of guest post on how to run AIMALL, which is our selected program for the integration algorithm.

The first thing we need is a WFN or WFX file, which contains the wavefunction in a Fortran unformatted file on which the Laplacian integration is to be performed. This is achieved in Gaussian09 by incluiding the keyword output=wfn or output=wfx in the route section and adding a name for this file at the bottom line of the input file, e.g.

filename.wfn

(NOTE: WFX is an eXtended version of  WFN; particularly necessary when using pseudopotentials or ECP’s)

Analyzing this file requires the use of a third party software such as AIMALL suite of programs, of which the standard version is free of charge upon registration to their website.

OpenAIMStudio (the accompanying graphical interface) and select the AIMQB program from the run menu as shown in figure 1.

 

Figure 1

Figure 1

Select your WFN/WFX file on which the calculation is to be run. (Figure 2)

 

Figure 2

Figure 2

You can control several options for the integration of the Laplacian of the electron density as well as other features. If your molecules are simple enough, you may go through with a successful and meaningful calculation using the default settings. After the calculation is finished, several result files are obtained. We’ll work in this tutorial only with *.mpgviz (which contains information about the molecular graph, MG) and *.sum (which contains all of  needed numerical data).

Visualization of the MG yields different kinds of critical points, such as: 1) Nuclear Attractor Critical Points (NACP); 2) Bond Critical Points (BCP); 3) Ring CP’s (RCP); and 4) Cage CP’s (CCP).

Of the above, BCP are the ones that indicate the presence of a chemical bond between two atoms, although this conclusion is not without controversy as pointed out by Foroutan-Njead in his paper: C. Foroutan-Nejad, S. Shahbazian and R. Marek, Chemistry – A European Journal, 2014, 20, 10140-10152. However, at a first approximation, BCP’s can help us to explore chemical interactions.

Now, let’s go back to visualizing those MGs (in our examples we’ve used methane and ethylene and acetylene). We open the corresponding *.mpgviz file in AIMStudio and export the image from the file menu and using the save as picture option (figure 3).

Figure 3

Figure 3

The labeled atoms are NACP’s while the green dots correspond to BCP’s. Multiplicity of a bond cannot be discerned within the MG; in order to find out whether a bond is a single, double or triple bond we have to look into the *.sum file, in which we’ll take a look at the bond orders between pairs of atoms in the section labeled “Diatomic Electron Pair Contributions and Delocalization Data” (Figure 4).

Figure 4

Figure 4

Delocalization indexes, DI’s, show the approximate number of electrons shared between two atoms. From the above examples we get the following DI(C,C) values: 1.93 for C2H4 and 2.87 for C2H2; on the other hand, DI(C,H) values are  0.98 for CH4, 0.97 in C2H4 and 0.96 in C2H2. These are our usual bond orders.

This is the first part of a crash tutorial on AIM, in my opinion this is the very basics anyone needs to get started with this interesting and widespread method. Thanks to all who asked about QTAIM, now you have your long answer.

Thanks a lot to my good friend Dr Pezhman Zarabadi-Poor for providing this contribution to the blog, we hope you all find it helpful. Please share and comment.

If a .fchk file wont open in GaussView5.0


I’ve found the following error regarding the opening of .fchk files in GaussView5.0.

CConnectionGFCHK::Parse_GFCHK()
Missing or bad data: Alpha Orbital Energies
Line Number 1234

The error is prevented to a first approximation (i.e. it at least will allow GV to open and visualize the file but other issues may arise) by opening the file and modifying the number of basis functions to equal the number of independent functions (which is lower)

FILE HEADER 
FOpt RM062X 6-311++G(d,p) 
Number of atoms I 75
Info1-9 I N= 9
 163 163 0 0 0 110
 2 18 -502
Charge I 0
Multiplicity I 1
Number of electrons I 314
Number of alpha electrons I 157
Number of beta electrons I 157
Number of basis functions I 1199
Number of independent functions I 1199
Number of point charges in /Mol/ I 0
Number of translation vectors I 0
Atomic numbers I N= 75
... ...
... ...

Once both numbers match you can open the file normally and work with it. My guess is this will continue to happen with highly polarized basis sets but I need to run some tests.

Transition State Search (QST2 & QST3) and IRC with Gaussian09


Theoretical evaluation of a reaction mechanism is all about finding the right transition states (TS) but there are no guarantees within the available methods to actually find the one we need. Chemical intuition in the proposal of a mechanism is paramount. Let’s remember that a TS is a critical point on a Potential Energy Surface (PES) that is a minimum in every dimension but one. For a PES with more than two degrees of freedom, a hyper-surface, envisioning the location of a TS is a bit tricky, in the case of a three dimensional PES (two degrees of freedom) the saddle point constitutes the location of the TS as depicted in figure 1 by a section of a revolution hyperboloid.

400px-Saddle_point

Fig1. Saddle point on a surface (min in one direction; max in the other)

Fig 1a Pringles chips -Yuck-. They exhibit a maximum on the direction parallel to the screen and a minimum on the direction perpendicular to the screen at the same point.

Fig 1a Pringles chips -Yuck-. They exhibit a maximum on the direction parallel to the screen and a minimum on the direction perpendicular to the screen at the same point.

The following procedure considers gas phase calculations. Nevertheless, the use of the SCRF keyword activates the implicit solvent calculation of choice in order to evaluate to some degree the solvent influence on the reaction energetics at different temperatures with the use of the temperature keyword.

The first step consists of a high level optimization of all minimums involved, such as reagents, products and intermediates, with a subsequent frequency analysis that includes no imaginary eigenvalues.

In order to find the structures of the transition states we use in Gaussian the Synchronous Transit-guided Quasi-Newton method [1] through the keywords QST2 or QST3. In the former case, coordinates for the reagents and products are needed as input; for the latter keyword, coordinates for the TS structure guess is needed also.

QST2)

%chk=file.chk
%nprocshared=n
%mem=nGB

#p opt=(qst2,redundant) m062x/6-31++G(d,p) freq=noraman Temperature=373.15 SCRF=(Solvent=Water)

Title card for reagents

Q M
Cartesian Coordinates for reagents
blank line
Title card for products

Q M
Cartesian Coordinates for products
blank line

QST3)

%chk=file.chk
%nprocshared=n
%mem=nGB

#p opt=(qst3,redundant) m062x/6-31++G(d,p) freq=noraman Temperature=373.15 SCRF=(Solvent=Water)

Title Card for reagents

Q M
Cartesian Coordinates for reagents
blank line
Title card for products

Q M
Cartesian Coordinates for products
blank line–
Title card for TS
Q M
Cartesian Coordinates for TS
blank line

NOTE: It is fundamental that the numbering order is kept constant throughout the molecular specifications of all two, or three, input structures. Hence, I recommend to build a set of molecules, save their structure, and then modified the coordinates on the same file to produce the following structure, that way the number for every atom will remain the same for every step.

As I wrote above, there are no guarantees of finding the right TS so many attempts are probably needed. Once we have the optimized structures for all the species involved in our mechanistic proposal we can plot their energies very simply with MS Excel the way we’ve previously described in this blog (reblogged from eutactic.wordpress.com)

Once we’ve succeeded in finding the structure of our TS we may run an Internal Reaction Coordinate (IRC) calculation. This calculation will connect the TS structure to those of the products and the reagents. Initial constant forces are required and these are commonly retrieved from the TS calculation checkpoint file through the RCFC keyword.

%chk=QST3_2p.chk
%nprocshared=8

#p m062x/6-31++G(d,p) IRC=(Maxpoints=50,RCFC,phase=(2,1))Temperature=373.15 SCRF=(Solvent=Water) geom=allcheck

Title Card

Q M
blank line

Finally, the IRC path can be visualized with GaussView from the Results menu. A successful IRC will link both structures along a single reaction coordinate proving that both reagents and products are linked by the obtained TS.

Hat tip to Howard Diaz who has become quite skillful in calculating these mechanisms as proven by his recent poster at the XII RMFQT a couple of weeks back. And as usual thanks to everyone who reads, comments, likes, recommends, rates and shares my silly little posts.

XIIth Mexican Reunion on Theoretical Physical Chemistry


As every year this month we had the yearly Mexican Reunion on Theoretical Physical Chemistry organized by prominent researchers in the field, such as Dr. Emilio Orgaz (UNAM), Dr. Alberto Vela (CINVESTAV) and many other. Over 150 different works were presented during this edition which took place in Juriquilla, Querétaro at one of the many campuses of the National Autonomous University of Mexico scattered all around the country. Below you can see some pictures from the talks and the first poster session.

20131119-154119.jpg20131119-154102.jpg20131119-154222.jpg

This time we contributed with a small poster on a mechanism proposed by Howard Diaz (an undergrad student from UAEM) on the equilibrium transformation of dihydrocinolines into 1-amino-indoles by an intramolecular rearrangement. May this post also serve as the starting point of a -mini-tutorial on how to evaluate a mechanism theoretically using QST3 and IRC in implicitly solvated environments (PCM)

20131119-154028.jpg

Howard Diaz posing next to his poster

The equilibrium under study and the proposed mechanism  by which it occurs, originally proposed by Frontana-Uribe et al. looks a bit like this:

equilibrium

Dihydrocinolines in equilibrium with 1-aminoindole

mechanism

Mechanistic proposal by Frontana-Uribe et al.

The energy profile, in which all transition states were calculated with the QST3 method, is presented below, calculated at various levels of theory. Also, the Internal Reaction Coordinate (IRC) connecting both states was calculated and is shown further below in the full poster.

Energy Profile

Energy Profile

From this results we believe that a new mechanistic proposal is needed since the energy barrier for the first step is quite high (~60 kcal/mol) and hence a bit unlikely to occur through that transition state. Nevertheless this is a first approach to elucidating a mechanism and the more knowledge about it the higher the control will be on this chemical transformation.

A full version of the poster is shown below for your convenience (Spanish). See you all at the next RMFQT in Morelia 2014!

Full Poster

Full Poster

Follow

Get every new post delivered to your Inbox.

Join 1,810 other followers

%d bloggers like this: