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.
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  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.
#p opt=(qst2,redundant) m062x/6-31++G(d,p) freq=noraman Temperature=373.15 SCRF=(Solvent=Water)
Title card for reagents
Cartesian Coordinates for reagents
Title card for products
Cartesian Coordinates for products
#p opt=(qst3,redundant) m062x/6-31++G(d,p) freq=noraman Temperature=373.15 SCRF=(Solvent=Water)
Title Card for reagents
Cartesian Coordinates for reagents
Title card for products
Cartesian Coordinates for products
Title card for TS
Cartesian Coordinates for TS
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.
#p m062x/6-31++G(d,p) IRC=(Maxpoints=50,RCFC,phase=(2,1))Temperature=373.15 SCRF=(Solvent=Water) geom=allcheck
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.
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.
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)
The equilibrium under study and the proposed mechanism by which it occurs, originally proposed by Frontana-Uribe et al. looks a bit like this:
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.
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!
It’s been a long time since I last posted something and so many things have happened in our research group! I should catch up with them in short but times have just been quite hectic.
Here is a contribution from Igor Marques at the University of Aveiro in Portugal (Group Website); the original text can be found as a comment in the original NBO Visualization post but it is pretty much the same thing you can find in this post. Here is a link to Chemcraft’s website. Thanks for sharing this, Igor!
=> Examples provided by Igor Marques used Chemcraft Version 1.7, build 365 <=
In the Gaussian input, with the NBORead option included under the population keyword, we should include the PLOT option as illustrated below. The gfoldprint keyword will print the basis set to the output file in the old G03 format. Some visualization programs require a certain format of the basis set to be printed to the output file in order to plot orbitals and other surfaces like the electron density; therefore, if you want to play safe, use gfoldprint, gfprint and gfinput in the same line. gfprint will print the basis set as a list but in the new G09 format, whereas gfinput will print the basis set using Gaussian’s own input format. (The used level of theory and number of shared processors are shown as illustrations only; also the Opt keyword is not fundamental to the visualization of the NBO’s)
%chk=filename.chk %nprocshared=8 #P b3lyp/6-311++g** Opt pop=(full,nboread) gfoldprint filename 0 1 molecular coordinates $NBO BNDIDX PLOT $END
this will generate files from *.31 to *.41
For the visualization of NBOs, you’ll need FILE.31 and FILE.37. Open FILE.31 from chemcraft. It will automatically detect FILE.37 (if in the same directory).
Tools > Orbitals > Render molecular orbitals
select the NBOs of interest (whcih are in the same order of the output),
Adjust settings > OK
On the left side of the window, select the NBO of interest and then click on ‘show isosurface’. Adjust the remaining settings. To represent another orbital, click on ‘keep this surface’ and then select another orbital from the rendered set and follow the previous steps.
> It’s possible to open a formated checkpoint file, containing the NBOs, in chemcraft.
%Chk=filename.chk %nprocshared=4 #P b3lyp/6-311++g** Opt pop=(full,nboread,savenbo) gfoldprint filename 0 1 molecular coordinates $NBO BNDIDX $END
the procedure is identical, but it is only necessary to read the *fchk file and then render the desired orbitals.
However, two problems might arise:
a) Orbitals in the checkpoint are reordered, thus requiring some careful inspection of the output.
b) Sometimes, for a larger molecule, the checkpoint might not be properly saved and the Gaussian job (as previously reported – http://goo.gl/DrSgA ) will end with:
Failed in SchOr1 in NBStor.
Error termination via Lnk1e in /data/programs/g09/l607.exe at Wed Mar 6 15:27:33 2013.
As usual, thanks to all for reading/commenting/rating this and other posts in this blog!
How to calculate the Delta G of solvation? This is a question that I get a lot in this blog, so it is about time I wrote a (mini)post on it, and at the same time put an end to this posting drought which has lasted for quite a few months due to a lot of pending work with which I’ve had to catch up. Therefore, this is another post in the series of SCRF calculations that are so popular in this blog. For the other posts on this subjects remember to click here and here.
SMD is the keyword you want to use when performing a Self Consistent Reaction Field (SCRF) calculation with G09. This keyword was only made available in this last version of the program and it corresponds to Truhlar’s and coworkers solvation model which is recommended by Gaussian itself as the preferred model to calculate Delta G of solvation. The syntax used is the standard way used in any other Gaussian input files as follows:
# 'route section keywords' SCRF=SMD
Separately, we must either perform a gas phase calculation or use the DoVacuum keyword within the same SCRF input, and then take the energy difference between gas phase and solvated models.
# 'route section keywords' SCRF=(SMD,DoVacuum)
No solvation or cavity model should be defined since, by definition, SMD will use the IEFPCM model which is a synonym for PCM.
As opposed to the previous versions of Gaussian, the output energy already contains all corrections, this is why we must take the difference between both values (remember to calculate them both at the same level of theory if calculated separately!). Nevertheless, when using the SMD keyword we get a separate line, just below the energy, stating the SMD-CDS non electrostatic value in kCal/mol.
The radii were also defined in the original paper by Truhlar; I’m not sure if using the keyword RADII with any of its options yields a different result or if it even ends in an error. Its worth the try!
Some calculation variations are not available when using SMD, such as Dis (calculation of the solute-solvent dispersion interaction energy), Rep (solute-solvent repulsion interaction energy) and Cav (inclusion of the solute cavitation energy in the total energy). I guess the reason for this might be that the SMD model is highly parametrized.
Have you found any issue with any item listed above? Pleases share your thoughts in the comments section below. As usual I hope this post was useful and that you all rate it, like it and comment.
I worked out a very quick and easy way to generate level schemes in Excel, based on a query from one of the other students in the group. Normally I would resort to something like the astonishing TikZ for this sort of task, however our group is very much a Microsoft Office 'What You See Is A Metaphor For…
I am frequently asked how to include an extra set of basis functions in a calculation or how to use an entirely external basis set. Sometimes this question also implies the explicit declaration of an external pseudopotential or Effective Core Potential (ECP).
New basis sets and ECPs are published continuously in specialized journals all the time. The same happens with functionals for DFT calculations. The format in which they are published is free and usually only a list of coefficients and exponents are shown and one has to figure out how to introduce it in ones calculation. The EMSL Basis Set Exchange site helps you get it right! It has a clickable periodic table and a list of many (not all) different basis sets at the left side. Below the periodic table there is a menu from which one can select which program we want our basis set for; finally we click on “get basis set” and a pop-up window shows the result in the selected format along with the corresponding references for citation. A multiple query can be performed by selecting more than one element on the table, which generates a list that almost sure can be used as input without further manipulations. Dr. David Feller is to be thanked for leading the creation of this repository. More on the history and mission of the EMSL can be found on their About page. Because of my experience, the rest of the post addresses the inclusion of external basis sets in Gaussian, other programs such as NwChem will be addressed in a different post soon.
The correct format for inclusion of an external basis set is exemplified below with the inclusion of the 3-21G basis set for Carbon as obtained from the EMSL Basis Set Exchange site (blank lines are marked explicitly just to emphasize their location:
spin multiplicity Molecular coordinates - blank line - C 0 S 3 1.00 172.2560000 0.0617669 25.9109000 0.3587940 5.5333500 0.7007130 SP 2 1.00 3.6649800 -0.3958970 0.2364600 0.7705450 1.2158400 0.8606190 SP 1 1.00 0.1958570 1.0000000 1.0000000 **** - blank line -
The use of four stars ‘****’ is mandatory to indicate the end of the basis set specification for any given atom. If a basis set is to be declared for a second atom, it should be included after the **** line without any blank line in between.
WARNING! Sometimes we can find more than one basis set in a single file this is due to different representations, spherical or cartesian basis sets. Gaussian by default uses cartesian (5D,7F) functions. Pure gaussian use 6 functions for d-type orbitals and 10 for f-type orbitals (6D, 10F). Calculations must be consistent throughout, hence all basis functions should be either cartesian or pure.
Inclusion of a pseudopotential allows for more computational resources to be used for calculation of the electronic structure of the valence shell by replacing the inner electrons for a set of functions which simulate the presence of these and their effect (such as shielding) on the valence electrons. There are full core pseudopotentialas, which replace the entire core (kernel). There are also medium core pseudopotentials which only replace the previous kernel to the full one, allowing for the outermost core electrons to be explicitly calculated. The correct inclusion of a pseudopotential is shown below exemplified by the LANL2DZ ECP by Hay and Wadt for the Chlorine atom.
spin multiplicity Molecular coordinates - blank line - basis set for atom1 **** basis set for atom2 (if there is any) **** - blank line - CL 0 CL-ECP 2 10 d potential 5 1 94.8130000 -10.0000000 2 165.6440000 66.2729170 2 30.8317000 -28.9685950 2 10.5841000 -12.8663370 2 3.7704000 -1.7102170 s-d potential 5 0 128.8391000 3.0000000 1 120.3786000 12.8528510 2 63.5622000 275.6723980 2 18.0695000 115.6777120 2 3.8142000 35.0606090 p-d potential 6 0 216.5263000 5.0000000 1 46.5723000 7.4794860 2 147.4685000 613.0320000 2 48.9869000 280.8006850 2 13.2096000 107.8788240 2 3.1831000 15.3439560
If a second ECP is to be introduced, it should be placed right after the first one without any blank line! If a blank line is detected then the program will assume it’s done reading all ECPs and Basis Sets.
Finally, here is an example of a combination of both keywords. If a second ECP was needed then we’d place it at the end of the first one without a blank line. The molecule is any given chlorinated hydrocarbon (H, C and Cl atoms exclusively)
#P B3LYP/gen pseudo=read ADDITIONAL-KEYWORDS - blank line - 0 1 Molecular Coordinates - blank line - H 0 S 3 1.00 19.2384000 0.0328280 2.8987000 0.2312040 0.6535000 0.8172260 S 1 1.00 0.1776000 1.0000000 **** C 0 S 7 1.00 4233.0000000 0.0012200 634.9000000 0.0093420 146.1000000 0.0454520 42.5000000 0.1546570 14.1900000 0.3588660 5.1480000 0.4386320 1.9670000 0.1459180 S 2 1.00 5.1480000 -0.1683670 0.4962000 1.0600910 S 1 1.00 0.1533000 1.0000000 P 4 1.00 18.1600000 0.0185390 3.9860000 0.1154360 1.1430000 0.3861880 0.3594000 0.6401140 P 1 1.00 0.1146000 1.0000000 **** Cl 0 S 2 1.00 2.2310000 -0.4900589 0.4720000 1.2542684 S 1 1.00 0.1631000 1.0000000 P 2 1.00 6.2960000 -0.0635641 0.6333000 1.0141355 P 1 1.00 0.1819000 1.0000000 **** - blank line - CL 0 CL-ECP 2 10 d potential 5 1 94.8130000 -10.0000000 2 165.6440000 66.2729170 2 30.8317000 -28.9685950 2 10.5841000 -12.8663370 2 3.7704000 -1.7102170 s-d potential 5 0 128.8391000 3.0000000 1 120.3786000 12.8528510 2 63.5622000 275.6723980 2 18.0695000 115.6777120 2 3.8142000 35.0606090 p-d potential 6 0 216.5263000 5.0000000 1 46.5723000 7.4794860 2 147.4685000 613.0320000 2 48.9869000 280.8006850 2 13.2096000 107.8788240 2 3.1831000 15.3439560 - blank line -
If you like this post or found it useful please leave a comment, share it or just give it a like. It is as much fun to find out people is reading as it is finding the answer to ones questions in someone else’s blog
One of the most successful posts this blog has ever published was on certain nuances of the solvation calculations on PCM in G03. However there are some differences in the SCRF modules between G09 and G03 and I here present some of them as well as some tips to work with the new version.
The SCFVAC keyword used to calculate the Gibbs Solvation Energy change is no longer available. It is now replaced by DoVacuum which should be included in the SCRF options as SCRF=(DoVacuum,etc.). However, the absolute solvation energy now requires a gas-phase optimization along with a frequency calculation followed by the same calculations with the SCRF=SMD option in the desired solvent and with the appropriate variables.
Gaussian 03 used to calculate and report non-electrostatic contributions to the solvation energy, however they were not included in the total energy nor during optimization procedures. These non-electrosatic interactions are no longer calculated in the default. In order to include these terms during the SCF procedure, and to have them reported separately, the SCRF=SMD option should be used.
My previous post on PCM mentioned the usage of the options OFac=0.8 & RMin=0.5 as part of the additional input. These ‘magic numbers’ (I hate the term) were used to modify the way by which the overlapping spheres were treated in order to create the surface which in turn defined the cavity. G09 uses a new algorithm to make the overlaps generate a smoother surface. I recommend to use the default values before including ‘magic parameters’.
All the default values which G03 used can be retrieved with the G03Defaults keyword, but it is strongly suggested to use it only for comparison with calculations previously done with the older version.
As with some other so-called ‘white papers’ this post will be further improved as more information arises during my own calculations. Thanks for reading! Please comment/like/share this post, as well as others in the blog, if you found useful the information within.
The Natural Bond Orbitals Deletion analysis provides an excellent approach to the assessment of bonding energy within a single molecular fragment or between many. It deletes specific elements of the Fock matrix (this means it sets their values to 0.000) and then re-diagonalizes it in order to find the difference in energy respect to the original matrix. About nine different kinds of deletions are available, which will be briefly summarized in the following section.
One of the main strengths of the NBO derived methods is their almost complete basis set independence, which allows us to obtain comparable numbers under different levels of theory.
Both G03 and G09 use the NBO3.1 program. The 5.0 version is sold separately by their creators, namely prof. Frank Weinhold, who can be contacted through their website. It’s not available for geometry optimizations (gradients),
some people insist on trying to get a different geometry by eliminating a certain interaction and that is just not possible directly with the NBODel method It is indeed possible to perform NBODeletions along with optimizations in G09 (Thanks to prof. Weinhold for his clarifying message) but there are some restrictions: molecular coordinates should be in Z-Matrix format and the number of variables to be optimized should not exceed 50; prof. Weinhold also recommends to use print=0 in the $NBO keylist in order to prevent the output files to become too big. Be sure to start with a proper geometry (close to the desired minimum) since given the nature of this analysis some apretiable geometry effects are to be expected.
The general syntax for its usage includes the string pop=NBODel in the route section of the GaussianX input file. Then, at the end of the file, the following is required:
--End of Input File-- --blank line-- $NBO $END $DEL Interactions to be deleted $END
ENTIRE BLOCKS OF ATOMS
In this kind of deletion one is able to delete all the elements between specific groups of atoms, as if their orbitals (and hence their common Fock elements) did not overlap.
--End of Input File-- --blank line-- $NBO $END $DEL ZERO 2 ATOM BLOCKS 2 BY 3 1 2 3 4 5 3 BY 2 3 4 5 1 2 $END --blank line--
The first line after $DEL indicates how many groups of atoms will be set to zero and the following lines indicate how many atoms belong to each group (i.e. the size of each block which in this case are 2 and 3, respectively). After this line the groups of atoms are listed, in this example all elements from atoms 1 and 2 with those of atoms 3, 4 and 5 will become zero. The next three lines are used for symmetry, so all the interactions from (1,2)->(3,4,5) are deleted along with (3,4,5)->(1,2)
DELETIONS BETWEEN ENTIRE MOLECULAR FRAGMENTS (Intermolecular deletions)
If we want to assess the interaction energy between two molecules, the previous method would consume a lot of time in declaring the size of each block with every atom of each molecule in it, plus there seems to be a limit to the size of the block. In this kind of deletion one is able to delete all the elements between two or more molecular fragments.
--End of Input File-- --blank line-- $NBO $END $DEL ZERO 2 DELOC FROM 1 TO 2 FROM 2 TO 1 $END --blank line--
The delocalizations can also be calculated only in one direction (FROM 1 to 2), in the case above both interactions 1->2 and 2->1 have been deleted. The input for a trimer in which all three fragments interacted with each other would look like this:
ZERO 6 DELOC FROM 1 TO 2 FROM 2 TO 1 FROM 2 TO 3 FROM 3 TO 2 FROM 1 TO 3 FROM 3 TO 1
In short, the number of bilateral delocalizations to be deleted is equal to twice the number of edges in a graph depicting the intermolecular interactions (A post on topology in chemistry is now due).
Reading the output file
Almost at the very end of the output file the following section can be found:
>>>>>>>>>> Convergence criterion not met.
SCF Done: E(RHF) = -4728.57245403 A.U. after 2 cycles
Convg = 0.2354D-03 -V/T = 2.0012
Energy of deletion : -4728.572454034
Total SCF energy : -4728.604640956
Energy change : 0.032187 a.u., 20.198 kcal/mol
The warning about the convergence can be disregarded without any concern about the accuracy of the outcome and it will show in every $DEL calculation. The SCF energy displayed in the second line is the energy corresponding to the modified Fock Matrix, which is the same as the one labeled as Energy of deletion. The Total SCF energy corresponds to the original Fock Matrix; the difference between them is labeled as Energy change and the value is reported in both atomic units as well as kcal/mol.
Some common errors and possible solutions
–> Sometimes you get the following error message at the beginning of the calculation making it crash:
** ERROR IN INITNF. NUMBER OF VARIABLES ( 57) **
** INCORRECT (SHOULD BE BETWEEN 1 AND 50) **
I have found that changing the molecule specification section from Z-matrix to Cartesian coordinates, or vice versa, overcomes this difficulty. Also, if the Opt keyword appears in the route section the previous message will be shown. Opt is
not available under the NBODel method (read the first paragraph for the proper correction).
–> Possible conflicts between NBODel and the usage of DFT methods:
In some revisions of Gaussian 03 there is a conflict when using NBODel and DFT methods. The IOp(5/48=10000) should be included to repair such conflict. This issue was solved in some revision of Gaussian 03 but I don’t know which, so try this if you have problems. Gaussian09 has taken care of the issue although still the usage of DFT to obtain NBODel calculations is not advised.
–> The following error is not self-explanatory:
Error termination via Lnk1e in ‘/../../path’
This particular error arises from the absence of the ‘$NBO $END’ line before the $DEL instruction. The previous line may or may not include additional keywords. If you are interested in computing some kind of deletion energy just leave the line as presented above in all previous examples. My guess is that the $DEL instruction does not calculate the corresponding NBO’s from which to make the deletion but it rather takes all the results from the $NBO instruction and works from there. Bottom line: don’t forget this line!
As with other posts tagged as ‘white papers’, this one will be updated and expanded every time new information is found. In the mean time, thanks to everyone for reading, commenting and rating, this keeps me going with the blog. Have you encountered problems with NBODel methods? share your experiences and solutions with the rest in the comments section.
Have a nice day!