# Category Archives: Gaussian

## Useful Thermochemistry from Gaussian Calculations

Statistical Mechanics is the bridge between microscopic calculations and thermodynamics of a particle ensemble. By means of calculating a partition function divided in electronic, rotational, translational and vibrational functions, one can calculate all thermodynamic functions required to fully characterize a chemical reaction. From these functions, the vibrational contribution, together with the electronic contribution, is the key element to getting thermodynamic functions.

Calculating the Free Energy change of any given reaction is a useful approach to asses their thermodynamic feasibility. A large negative change in Free Energy when going from reagents to products makes up for a quantitative spontaneous (and exothermic) reaction, nevertheless the rate of the reaction is a different story, one that can be calculated as well.

Using the **freq** option in your route section for a Gaussian calculation is mandatory to ascertain the current wave function corresponds to a minimum on a potential energy hypersurface, but also yields the thermochemistry and thermodynamic values for the current structure. However, thermochemistry calculations are not restricted to minima but it can also be applied to transition states, therefore yielding a full thermodynamic characterization of a reaction mechanism.

A regular **freq** calculation yields the following output (all values in atomic units):

Zero-point correction= 0.176113 (Hartree/Particle) Thermal correction to Energy= 0.193290 Thermal correction to Enthalpy= 0.194235 Thermal correction to Gibbs Free Energy= 0.125894 Sum of electronic and zero-point Energies= -750.901777 Sum of electronic and thermal Energies= -750.884600 Sum of electronic and thermal Enthalpies= -750.883656Sum of electronic and thermal Free Energies= -750.951996

For any given reaction say A+B -> C one could take the values from the last row (lets call it G) for all three components of the reaction and perform the arithmetic: DG = GC – [GA + GB], so products minus reagents.

By default, Gaussian calculates these values (from the previously mentioned partition function) using normal conditions, T = 298.15 K and P = 1 atm. For an assessment of the thermochemistry at other conditions you can include in your route section the corresponding keywords **Temperature=**x.x and **Pressure=**x.x, in Kelvin and atmospheres, respectively.

(Huge) **Disclaimer**: Although calculating the thermochemistry of any reaction by means of DFT calculations is a good (and potentially very useful) guide to chemical reactivity, getting quantitative results require of high accuracy methods like G3 or G4 methods, collectively known as Gn mehtods, which are composed of pre-defined stepwise calculations. The sequence of these calculations is carried out automatically; no basis set should be specified. Other high accuracy methods like CBS-QB3 or W1U can also be considered whenever Gn methods are too costly.

## Atom specifications unexpectedly found in input stream.

*“Well, where else were they supposed to appear?”*

I was sent this error along with the previous question for a failed optimization. Apparently there is no answer in the internet (I quickly checked) so here it is:

Gaussian is confused about finding atomic coordinates because there is also a** geom=check **instruction placed in the route section, i.e., it was told to retrieve the atomic coordinates from a checkpoint and then it was given those atomic coordinates within the input so it doesn’t know what you mean and exits.

## Calculating NMR shifts – Short and Long Ways

Nuclear Magnetic Resonance is a most powerful tool for elucidating the structure of diamagnetic compounds, which makes it practically universal for the study of organic chemistry, therefore the calculation of ^{1}H and ^{13}C chemical shifts, as well as coupling constants, is extremely helpful in the assignment of measured signals on a spectrum to an actual functional group.

Several packages offer an additive (group contribution) empirical approach to the calculation of chemical shifts (ChemDraw, Isis, ChemSketch, etc.) but they are usually only partially accurate for the simplest molecules and no insight is provided for the more interesting effects of long distance interactions (*vide infra*) so quantum mechanical calculations are really the way to go.

With Gaussian the calculation is fairly simple just use the NMR keyword in the route section in order to calculate the NMR shielding tensors for relevant nuclei. Bear in mind that an optimized structure with a large basis set is required in order to get the best results, also the use of an implicit solvation model goes a long way. The output displays the value of the total isotropic magnetic shielding for each nucleus in ppm (image taken from the Gaussian website):

Magnetic shielding (ppm): 1 C Isotropic = 57.7345 Anisotropy = 194.4092 XX= 48.4143 YX= .0000 ZX= .0000 XY= .0000 YY= -62.5514 ZY= .0000 XZ= .0000 YZ= .0000 ZZ= 187.3406 2 H Isotropic = 23.9397 Anisotropy = 5.2745 XX= 27.3287 YX= .0000 ZX= .0000 XY= .0000 YY= 24.0670 ZY= .0000 XZ= .0000 YZ= .0000 ZZ= 20.4233

Now, here is why this is the long way; in order for these values to be meaningful they need to be contrasted with a reference, which experimentally for ^{1}H and ^{13}C is tetramethylsilane, TMS. This means you have to perform the same calculation for TMS at -preferably- the same level of theory used for the sample and substract the corresponding values for either H or C accordingly. Only then the chemical shifts will read as something we can all remember from basic analytical chemistry class.

GaussView 6.0 provides a shortcut; open the Results menu, select NMR and in the new window there is a dropdown menu for selecting the nucleus and a second menu for selecting a reference. In the case of hydrogen the available references are TMS calculated with the HF and B3LYP methods. The SCF – GIAO plot will show the assignments to each atom, the integration simulation and a reference curve if desired.

The chemical shifts obtained this far will be a good approximation and will allow you to assign any peaks in any given spectrum but still not be completely accurate though. The reasons behind the numerical deviations from calculated and experimental values are many, from the chosen method to solvent interactions or basis set limitations, scaling factors are needed; that’s when you can ask the Cheshire Cat which way to go

If you don’t know where you are going any road will get you there.

Lewis Carroll – Alice in Wonderland

Well, not really. The Chemical Shift Repository for computed NMR scaling factors, with Coupling Constants Added Too (aka CHESHIRE CCAT) provides with straight directions on how to correct your computed NMR chemical shifts according to the level of theory without the need to calculate the NMR shielding tensor for the reference compound (usually TMS as pointed out earlier). In a nutshell, the group of Prof. Dean Tantillo (UC Davis) has collected a large number of isotropic magnetic shielding values and plotted them against experimental chemical shifts. Just go to their scaling factors page and check all their linear regressions and use the values that more closely approach to your needs, there are also all kinds of scripts and spreadsheets to make your job even easier. Of course, if you make use of their website **don’t forget** to give the proper credit by including these references in your paper.

We’ve recently published an interesting study in which the 1H – 19F coupling constants were calculated via the long way (I was just recently made aware of CHESHIRE CCAT by Dr. Jacinto Sandoval who knows all kinds of web resources for computational chemistry calculations) as well as their conformational dependence for some substituted 2-aza-carbazoles (fig. 1).

The paper is published in the Journal of Molecular Structure. In this study we used the GIAO NMR computations to assign the peaks on an otherwise cluttered spectrum in which the signals were overlapping due to conformational variations arising from the rotation of the C-C bond which re-orients the F atoms in the fluorophenyl grou from the H atom in the carbazole. After the calculations and the scans were made assigning the peaks became a straightforward task even without the use of scaling factors. We are now expanding these calculations to more complex systems and will contrast both methods in this space. Stay tuned.

## Error for Gaussian16 .log files and GaussView5

There’s an error message when opening some Gaussian16 output files in GaussView5 for which the message displayed is the following:

ConnectionGLOG::Parse_Gauss_Coord(). Failure reading oriented atomic coordinates. Line Number

We have shared some solutions to the GaussView handling of *chk and *.fchk files in teh past but never for *.log files, and this time Dr. Davor Šakić from the University of Zagreb in Croatia has brought to my attention a fix for this error. If “Dipole orientation” with subsequent orientation is removed, the file becomes again readable by GaussView5.

Here you can download a script to fix the file without any hassle. The usage from the command line is simply:

˜$ chmod 777 Fg16TOgv5 ˜$ ./Fg16TOgv5 name.log

The first line is to change and grant all permissions to the script (use at your discretion/own risk), which in turn will take the output file **name.log** and yield two more files: **gv5_name.log** and and **name.arch**; the latter archive allows for easy generation of SI files while the former is formatted for GaussView5.x.

Thanks to Dr. Šakić for his script and insight, we hope you find it useful and if indeed you do please credit him whenever its due, also, if you find this or other posts in the blog useful, please let us know by sharing, staring and commenting in all of them, your feedback is incredibly helpful in justifying to my bosses the time I spent curating this blog.

Thanks for reading.

## Python scripts for calculating Fukui Indexes

One of the most popular posts in this blog has to do with calculating Fukui indexes, however, when dealing with a large number of molecules, our described methodology can become cumbersome since it requires to manually extract the population analysis from two or three different output files and then performing the arithmetic on them separately with a spreadsheet or something.

Our new team member Ricardo Loaiza has written a python script that takes the three aforementioned files and yields a .csv file with the calculated Fukui indexes, and it even points out which of the atoms exhibit the largest values so if you have a large molecule you don’t have to manually check for them. We have also a batch version which takes all the files in any given directory and performs the Fukui calculations for each, provided it can find file triads with the naming requirements described below.

Output files must be named *filename*.log (the N electrons reference state), *filename***_plus**.log (the state with N+1 electrons) and *filename***_minus**.log (the N-1 electrons state). Another restriction is that so far these scripts only work with NBO population analysis as provided by the NBO3.1 program available in the various versions of Gaussian. I imagine the listing is similar in NBO5.x and NBO6.x and so it should work if you do the population analysis with them.

The syntax for the single molecule version is:

python fukui.py filename.log filename_minus.log filename_plus.log

For the batch version is:

./fukuiPorLote.sh

(*Por Lote* means *In Batch* in Spanish.)

These scripts are available via GitHub. We hope you find them useful, and you do please let us know whether here at the comments section or at our GitHub site.

## fchk file errors (Gaussian) Missing or bad Data: RBond

We’ve covered some common errors when dealing with formatted checkpoint files (*.fchk) generated from Gaussian, specially when analyzed with the associated GaussView program. (see here and here for previous posts on the matter.)

Prof. Neal Zondlo from the University of Delaware kindly shared this solution with us when the following message shows up:

CConnectionGFCHK::Parse_GFCHK() Missing or bad data: Rbond Line Number 1234

The Rbond label has to do with the connectivity displayed by the visualizer and can be overridden by close examination of the input file. In the example provided by Prof. Zondlo he found the following line in the connectivity matrix of the input file:

2 9 0.0

which indicates a zero bond order between atoms 2 and 9, possibly due to their proximity. He changed the line to simply

2

So editing the connectivity of your atoms in the input can help preventing the Rbond message.

I hope this helps someone else.

## Quick note on WFN(X) files and MP2 calculations #G09 #CompChem

A few weeks back we wrote about using WFN(X) files with MultiWFN in order to find σ-holes in halogen atoms by calculating the maximum potential on a given surface. We later found out that using a chk file to generate a wfn(x) file using the guess=(read,only) keyword didn’t retrieve the MP2 wavefunction but rather the HF wavefunction! Luckily we realized this problem very quickly and were able to fix it. We tried to generate the wfn(x) file with the following keywords at the route section

#p guess=(read,only) density=current

but we kept retrieving the HF values, which we noticed by running the corresponding HF calculation and noticing that every value extracted from the WFN file was exactly the same.

So, if you want a WFN(X) file for post processing an MP2 (or any other post-HartreFock calculation for that matter) ask for it from the beginning of your calculation in the same job. I still don’t know how to work around this or but will be happy to report it whenever I do.

PS. A sincere apology to all subscribers for getting a notification to this post when it wasn’t still finished.

## 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 MultiplicityCoordinates of reagents --blank line—Charge MultiplicityCoordinates of products --blank line--Charge MultiplicityCoordinates 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 testgeom=connectivity --blank line--Charge MultiplicityZ-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=calcfcgeom=connectivity --blank line--Charge MultiplicityCoordinates 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 MultiplicityCoordinates 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 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=qst2geom=connectivity freq=noraman --blank line--Charge MultiplicityCoordinates of reagents --blank line--Charge MultiplicityCoordinates 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

*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.*

**Connection.**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=tsgeom=connectivity freq=noraman --blank line--Charge MultiplicityCoordinates 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 MultiplicityCoordinates 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!