Blog Archives

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 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!

 

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.

Restarting Calculations from rwf Files – Gaussian


Having a long calculation terminated just because it ran out of time in the queue is very frustrating; even more so if restarting it from the last accesible point is hard to do.

I have recently performed some particularly demanding calculation: Basis Set Superposition Error (BSSE) with the Counterpoise method and second order Moller-Plesset perturbation theory calculation (MP2). The calculation ran out of time but I was able to restart it because I had the rwf file! My input looked a bit like this:

#p mp2/GEN counterpoise=2 maxdisk=200GB

So here is how it works.

The very first line of your calculation gives you the process ID number which is not necessarily the same as the PID given by the queue system (in fact, is not the same because the latter corresponds to the submitted script, not the instructions in it i.e. your calculation)

 Entering Gaussian System, Link 0=g09
 Initial command:
 /opt/SC/aplicaciones/g09-C.01/l1.exe /tmpu/joaqbf_g/joaqbf/Gau-38954.inp 
-scrdir=/tmpu/joaqbf_g/joaqbf/
 Entering Link 1 = /opt/SC/aplicaciones/g09-C.01/l1.exe PID=     38955.

(emphasis in red is mine)

This is the number you want to write down. You will need to find the corresponding rwf file (usually in your SCRATCH directory) as Gau-PID.rwf (in the aforementioned case, Gau-38955.rwf). If you are a bit paranoid like myself you want to copy and keep this file safe but be aware that these are very long files, in my case it was 175 GB long. Now you need to launch your calculation again with the following input:

%rwf=myfile.rwf
%nosave
%chk=myfile.chk

Title Card

# restart

rest of input

You can add all other controls to the Link0 section such as %nprocshared or %mem according to your needs.

I’m pretty sure it should work for other kinds of calculations in which taking from the checkpoint file is not as easy, so if you run into this kind of problems, its worth the try.

Some .fchk files wont open in GaussView5.0 (Update)


A couple of weeks ago I posted a solution for a common error regarding .fchk files that will display the error below when opened with GaussView5.0. As I expected, this error has to do with the use of diffuse functions in the basis set and is related to a change of format between Gaussian versions.

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

Although the method described in the previous post works just fine, the following update is a better approach. Due to a change of spelling between G03 and G09 (which has been corrected for G09 but not available for GV versions prior to 5.0.9) one must change “independent” for “independant

To make the change directly from the terminal the following command is needed:

sed -i 's/independent/independant/g' file.fchk

Alternatively you can redirect the output to a new file

sed -e 's/independent/independant/g' file.fchk > newfile.fchk

if you want to keep the old version and work with a new one.

Of course this edition can be performed manually with any text editor available (for example if you work in Windows) but solutions from the terminal always seem easier and a lot more fun to me.

Thanks to Dr. Fernando Cortés for sharing his insight into this issue.

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.


This is the first time I reblog a post from a fellow computational chemist and the reason why I do it is because of its beautiful simplicity and usefulness. Given the scope this blog has taken I think this post becomes most appropriate. This post will show you how to create an energy level diagram using nothing but MS Excel.
Kudos to ‘Eutactic’, from Australia, for coming up with a nice solution to this problem. Check out his blog at eutactic.wordpress.com.
Thanks for letting me repost it 🙂

eutactic

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 Cosmic Horror‘ group and recommending that a colleague learns two new markup languages to produce a figure is probably not helpful in the short term. One of the issues with charting energy levels in Excel is that levels are typically represented by horizontal bars connected at their vertices with lines representing transitions. Whilst Excel does have a horizontal bar as a marker, it possesses two show-stopping limitations:

  1. It is only uniformly scalable, and can only be scaled so far – we cannot make it anywhere near wide and…

View original post 222 more words

Some software trickery…


This post will become updated continuously in order to include as much of these useful tips as I find along the way, so if you are interested please subscribe so you don’t have to visit often. Some are old or even a bit obsolete in terms of software versions but I still include them so they can all be gathered in one place; this is basically a summary of some common errors in the use of Gaussian (G03 and G09). Please feel free to comment on this post, sharing your tricks with various programs in whatever branch of computational chemistry you might be working on.

Vibrations calculated with Gaussian 09 can’t be visualized with GaussView 3.x

*Change the lines “Atom  AN” with “Atom AN”. They look the same but there are TWO spaces in the first (this is the way G09 prints the output files) and in the second one there is only ONE space. The credit is entirely due to Jean Poully who posted this trick on the CCL a few months ago.

*Reader John Keller from Alaska suggests to use the formatted chk file (.fchk) to avoid this problem. As usual with formatted chk files make sure you format them in the same computer where you calculated it, or at least make sure both versions are compiled under the same architecture (32 or 64 bits).

Molecules appear flat on Molekel versions previous to 5.x

Molekel 4.x reads Gaussian98 files only. Gaussian03 and Gaussian09 files need to have the header corrected in the following part, right at the beginning of the file:

******************************************
Gaussian 09:  EM64L-G09RevA.02 11-Jun-2009
9-May-2010
******************************************

Change to:

******************************************
Gaussian 98:  EM64L-G09RevA.02 11-Jun-2009
9-May-2010
******************************************

A rather old one and kind of a weird one too but it works and some of us (myself very much included) like the previous Molekel far more than the new 5.x version.

Electrostatic potentials can not be visualized with Molekel 4.x

This is an error that stems from the fact that atomic charges are not read by Molekel, so all electrostatic properties derived from them are not calculated nor visualized. The problem is that Gaussian also changed the header for the charges when going from Gaussian98 to Gaussian03. In Gaussian 98 charges are labeled as “Total Atomic Charges”. In G03 and G09, atomic charges are labeled “Mulliken Atomic Charges”. Therefore find and change this label for the latter. PLEASE NOTE that Molekel will look for the title “Total Atomic Charges” so if your calculation was a geometry optimization then you have to find the set of charges that correspond to the last step of it. Look for the last string and change that one; you may also change them all with the right edition option in your editor, Molekel will use the last one since that should be the one corresponding to the optimized structure.

jbf

I’ll post something else when I find/remember it. Thanks for reading!

Mi amor!!!
Muchísimas felicidades!!! Estoy muy orgulloso de ti, de ser tu novio y tenerte en mi vida. Como persona y como profesional no eres más que admirable. Al igual que quienes te escriben, estoy seguro que este es solo el principio, mi amor, ese negocio en tus manos va a llegar hasta donde tu papá siempre quiso y aun más lejos. Te amo y estoy
%d bloggers like this: