F- + CH3Cl -> Cl- + CH3F
The general approach is as follows
1. Build the structure (F- + CH3Cl)
2. Optimize the geometry while constraining a key distance (between the F and C atom) to a certain value (2.0 Å)
3. Compute the frequencies for the optimized geometry to check for an imaginary frequency
4. If there is an imaginary frequency, use this geometry and frequency information as a starting point for the TS search
The distance constraint is imposed using
$zmat ifzmat(1)=1,1,6 fvalue(1)=2.0 $end
which is used in conjunction with the delocalized internal coordinates I discussed previously. 1,1,6 defines the distance (denoted with the first "1") between atoms 1 and 6. Note that this distance does not have be 2.0 Å in the guess geometry (though the closer it is to 2.0 Å the faster the optimization will converge).
I also show how the frequency information contained in the $HESS group, found at the end of the .dat file that GAMESS produces along with the output file, can be transferred to the input file for the TS search, where hess=read is used in $statpt. Any text editor can be used for this. Note that the hssend=.t. keyword has to be used to cause the frequencies to be computed at the end of the optimization. I have showed the use of this keyword in a screencast in a previous post.
The complete $statpt groups looks like this
$statpt opttol=0.0005 nstep=50 hess=read ihrep=5 hssend=.t. $end
ihrep=5 causes the frequencies to be recomputed at every 5th step. This really helps convergence but is computationally expensive.
Note that there are no constraints imposed when searching for the TS.
The key question is of course why I chose the C-F distance and 2.0 Å for the constraint value. This is simply the first thing that popped into my head and it happened to work. If it hadn't, I would have tried some other values and, if none worked, I would have tried to the C-Cl distance. Future posts will show examples of trying difference distances.
Because you may have to try different constraints before succeeding, it is a good idea to find the TS at a cheap level of theory, such as PM3. The C-F distance in the PM3 TS structure (2.10 Å) should provide a good guess for TS searches at higher levels of theory.
I tried it for this reaction using RHF/3-21G, but with no luck. That is because I forgot the very first step in searching for a TS: make sure there is one! In the screencast below I show how a RHF/3-21G geometry optimization (with no constraints) starting from F- + CH3Cl leads to Cl- + CH3F, meaning that the former structure is not a minimum on the gas phase RHF/3-21G potential energy surface. So there is not much point in looking for the corresponding TS.
Posts
35 comments:
Hi, your blogs nice.I'm studying chemical engineering and interested to gain skills of Molecular Modeling.Could I be able to learn them side by side doing my b.tech? And can you give some advices regarding the same.
Jan,
Thanks for the example! I've been busy running calcs and learning GAMESS and realized that I'm not confident in running conformer searches in GAMESS. Can you do a screencast that shows how to search for conformers using GAMESS? Thanks again.
NUchem
Begujam -
if you're completely new to molecular modeling, then my best advice is to find a course on it at your school. There are also some good books out there: Frank Jensen's "Introduction to Computational Chemistry", Chris Cramer's "Essentials of Computational Chemistry", and Andrew Leach's "Molecular Modeling: Principles and Applications". My book won't be out until next year ...
NUchem -
GAMESS (as far as I know) has no automatic conformer search option. You can use the conformer search option in Avogadro, and then create GAMESS geometry optimization input files for each unique geometry. Is that the kind of thing you had in mind?
Hi , you blog is amazing!
Now, I just working on a cycloaddition reaction, but I have some problems to find the transition state!!
The reaction is a Diels Alder, but I'm not understand in which way I have to constraint the bonds.
Can you help me please?
Best regards
Paolo
If we consider the reaction ethene + 1,3-butadiene -> cyclohexe, I suggest:
1) optimize cyclohexene
2) starting with this optimized structure, constrain the 2 C-C single bonds next to the double bond to 2 Å.
See if that gives you an imaginary frequency. If you haven't already done so, you should also look at the other posts labeled transition states.
Glad to hear you like the blog.
I saw the posts for TS, but are little bit different!!
In the point 2) i've to consider the cyclohexene alone or whit butadiene?
Sorry for the question!
I've made a stupid question, sorry :)
Sorry for disturbing.
I managed to find the transition state of simple reaction ethene + 1,3-butadiene -> cyclohexene;
now have to work on a more complex reaction: the product of the Diels-Alder has five six membered cycles condensed, one created by the reaction.
I tried to work by simplifying the molecule, removing the cycles non involved in the mechanism, and I managed to find the transition state.
However, when I work on the whole product, which impose constraints (the same as the simplified reaction) don't break the bonds.
The constraints are the same as the reaction ethene + 1,3-butadiene: in cyclohexene I impose the length of the two bonds that are formed in the reaction to 2.25A.
In which way I can obtain the TS?
Paolo
I am not sure what you mean when you say that the constraints don't break the bonds. Can you be more specific?
Two general pieces of advice:
Start with the TS you found without the rings, add the rings and minimize, while keeping the TS region frozen. This should be a better TS-guess than constraining the distances. This post shows how to do a constrained minimization in Avogadro.
Also, rather than adding all rings at once, start adding them one at a time, and use the TS you find for the smaller model as an initial guess for the next, larger, model.
Hope this helps.
When I say that the constraints don't break bonds, i means that product remains intact and doesn't give the imaginary vibrations. I think that is for the rigidity of the structure with all rings.
Now I try with your advice.
Thanks a lot for your willingness
I have add all the rings freezing the the TS region, but now in which way I have to use the new structure? I have to optimize the geometry or find another transition state?
Paolo
First check if this structure has an imaginary frequency with a mode that looks like what you want.
If it does, use it for a runtyp=sadpoint run.
Hi, I've been following your blog. I would like to learn about Molecular Modeling for my thesis. So I was assigned to make an IRC calculation on the reaction of Carnosine with 4-hydroxynonenal. Can help me out? What do I do first? Thanks in advance.
Ralph
By the way, here's the reaction that I will do. Thanks!
[IMG]http://i51.tinypic.com/2e2gx3q.jpg[/IMG]
Ralph
I have made a post on computing an IRC, which should hopefully help.
Hi Jan,
I have a question for you that makes me mad: I computed a transition state of a complex system; I got the imaginary frequency and the vectors seem concord to the mechanism. Now, plotting the Potential Energy Surface the calculated TS gives to me a negative activation energy. So the reaction seems to be barrierless. Is it possible? I know that usually increasing the temperature in real systems we can have a barrierless reaction but in a simulation seems strange. I used a PM3 method and I think that this probably is the reason. Performing the single point energy with a B3LYP/6311-G-(2d,p) the energy gap decrease but the TS energy remains 9 kcal lower than the reactants.
Help me please and thanks for your time
Paolo
If you found a proper TS, then its energy has to be higher than the reactant and product it connects (this does not account for single point energies).
But without knowing more it's hard to tell what's wrong:
1. What's the value of the imaginary frequency?
2. Do you use NVIB=2 in the Hessian calculation?
3. What's the value of the (negative) activation energy at the PM3 level?
4. Are you sure you found the lowest energy conformation of your reactant.
5. Have you computed an IRC?
Sure:
1. the value of the imaginary frequency is 225.3 cm^-1
2.yes
3. around 5 kcal/mol lower than reactant
4. yes
The pattern show a not barrierless reaction
link:http: //db.tt/dOfu7iAz
How I can explain this?
Paolo
I dont understand 4. Can you explain exactly how you obtained this plot?
I follow exactly the procedure that has been explained in your post concerning the Intrinsic Reaction Coordinate.
If thats the IRC why do you say "plotting the Potential Energy Surface the calculated TS gives to me a negative activation energy"?
To construct the Potential Energy Surface I use to optimize the reactants singly, the TS, and the product. After that I plot the three steps: in this way I get a negative activation energy and the reaction seems to be barrierless.
My idea is that in this case the sum of the optimized reactants,each one considered at infinity distance from the other one, is higher in energy respect to the TS.
I hope you understand now.....thanks for your help, really!!
Paolo
Sorry, I don't understand
:)
My way to build the PES:
Reactants---->TS---->Product
Reactants = react1 + react2
In your way to build IRC we use only the TS geometry so the reactants are considered interacting. In my way no. Is the simple sum of the energies.
However thanks a lot Prof.
P.S
After reading your post one year ago I learned to use GAMESS, I graduated myself with a magnificent Mass spectrometry thesis with some computational features (all learned by you!!!!!!)
Now I got a Research bourse in computational chemistry!!!
So thanks a lot!!
Now, I understand. You have
R1 + R2 -> R1/R2 -> TS -> P
and E(R1)+E(R2) > E(TS)
This can happen and there is nothing (necessarily) wrong with your calculations.
If you are trying to model a reaction in solution, then your activation energy is
Ea = E(TS)-E(R1/R2)
This is because the energy liberated by binding is lost to the environment before it can be used to react.
If you are trying to model a reaction in the gas phase, then you would say that this is a barrier-less reaction and the reaction rate is proportional to the collision frequency. This is because the energy cannot escape the complex, and can be used to overcome the barrier. However, this will be pressure dependent (higher pressure, worse assumption).
PS. I am really happy to hear you found the blog useful. Thanks for telling me.
:)
The discussion with Paolo lead to this post.
Thanks for these instructions. I used them to model a simple Sn2 reaction for an o chem class that I teach. I ran into some problems trying to model the reactants, and wanted to share the solution so that if someone else runs into the same problem, they know what to do.
When the two reactants were close together (<6Å or so), I could do minimization and energy computations with DFT just fine. But when the reactants were farther apart, turning DFT on made GAMESS abort because it could not get the SCF to converge. I eventually figured out that the problem was my DFT method. I was using B3LYP, which isn't designed for long-range correlations, like the separation between two reactants before their orbitals start to interact. Changing the DFT method to BLYP and turning on the long-range correlation correction in GAMESS by adding $DFT LC=.t. to the input file fixed the problem. Now the SCF converges just fine.
Hi,
I'm trying to copy the input that you wrote in GAMESS for optimizing the geometry for the reaction:
F- + CH3Cl -> Cl- + CH3F
The following error message always
results in the output: "ERROR: UNABLE TO PROJECT DLC".
Below is the input code that I've been using. What is the reason why I recieved the error message?
$CONTRL SCFTYP=RHF RUNTYP=OPTIMIZE
ICHARG=-1 MULT=1 COORD=ZMTMPC $END
$statpt opttol=0.0005 nstep=50 hssend=.t. $end
$force nvib=2 $end
$contrl nzvar=1 $end
$zmat dlc=.t. auto=.t. nonvdw(1)=1,6 1,2 $end
$zmat ifzmat(1)=1,1,6 fvalue(1)=2.0 $end
$BASIS GBASIS=PM3 $END
$DATA
CH3FCl(-1)
C1 1
C 0.0000000 0 0.0000000 0 0.0000000 0 0 0 0
Cl 1.7600005 1 0.0000000 0 0.0000000 0 1 0 0
H 1.0900004 1 109.47118 1 0.0000000 0 1 2 0
H 1.0899995 1 109.47124 1 -119.99998 1 1 2 3
H 1.0899995 1 109.47124 1 119.99998 1 1 2 3
F 2.0000001 1 171.80211 1 0.0000000 1 1 2 3
$END
In regards to the previous post, the input code that I typed into GAMESS did have a space before each $. The spaces did not appear when I copied the input into the blog post.
That's strange. It works fine for me, i.e. I copied the text from your comment, made an input file and ran it.
I suggest posting your question to the Google GAMESS group with lots of detail including version, hardware, and quotes from the output.
In orer to solve the error: UNABLE TO PROJECT DLC, I added IFDMOD=2 to the $zmat group. IFDMOD selects a method to freeze internals in DLC.
I also had to add AUTOFV=.f. into the $zmat group in order to allow GAMESS to read the fvalue input of 2.0. AUTOFV=.t. is the default and it tells GAMESS to generate an fvalue array automatically, thus ignoring values inputted by the user.
I have an outdated version of GAMESS. In the newer versions of GAMESS, IFDMOD is put in by default. I'm not certain, but I think that AUTO=.f. is also added a default in the newer versions.
OK, thanks for the follow up. Just curious, what version do you use. I use a fairly ancient version 12 JAN 2009 (R3) on my laptop.
I have PC GAMESS version 7.0 (14 August 2006).
Hello sir, I am trying to calculate TS energy of OH- + CH3Cl using your method. First I have not been able to get a negative frequency after several trials. Second, I don't know how you got the data about $HESS that you copied into your ts.inp file below the $DATA. Thanks in anticipation
Post a Comment