ChemWiki - User contributions [en]

Visualizing High Model Orbitals

2014-03-20T22:30:32Z

Lmt09: /* Aim */

==Visualizing High Model Orbitals==

===Aim===
When using ONIOM we are often interested in the effect of an environment on the model compound. For this reason it is often useful to visualize the orbitals of the high model region after a calculation. When trying to do this in Gaussview, however, it is found that the orbitals displayed are those of the low model. This is the reason why it was necessary to construct orbitals using ''guess=input'' in the ONIOM(CASSCF:AM1) examples. This tutorial explains how to access the orbitals after a calculation has been run, for example if we wish to localize orbitals to ensure the correct active space has been chosen in the previous example. In the case of QM:MM calculations it is possible to skip the punch orbitals step as there is only one set of orbitals on the checkpoint file; just get the model geometry and read the guess from the checkpoint file directly (I think newer versions of Gaussian actually die now if you try to do this - Lee).

===System===
In this tutorial we examine the spurious transition state in the diels-alder cycloaddition between maleic anhydride and cyclohexadiene. We extract the high model orbitals and localize them to ensure that the correct active space has been chosen. This allows us to check that the choice of active space is not the cause of disagreement with the high real calculation, which indicates a symmetric transition state.

[[Image:TS_male_cyc.jpg|frame|Spurious transition state of maleic anhydride and cyclohexadiene]]

===Method===
====Punch Orbitals====
The first task is to obtain the orbitals if the high model in a format that can be read back in by Gaussian. This can be achieved using the ''punch=MO'' keyword but, in order to punch the high model orbitals we need to use a nonstandard route.

First we use ''testrt'' to obtain the standard route:
<pre>
-----------------------------------------------------------------
#p oniom(casscf(6,6)/sto-3g:hf/sto-3g) guess=read nosymm punch=MO
-----------------------------------------------------------------
1/38=1,52=2/1;
2/12=2,15=1,17=6,18=5,40=1/2;
1/38=1,52=2,53=3172/20;
3/6=3,11=9,16=1,25=1,30=1,116=-2/1,2,3;
4/5=1,17=6,18=6/1;
5/5=2,38=6/2;
6/7=2,8=2,9=2,10=2,28=1/1;
1/52=2,53=2032/20;
3/6=3,16=1,25=1,32=1,116=101/1,2,3;
4/5=1,17=6,18=6/1,5;
5/5=2,17=1000000,38=6/10;
6/7=2,8=2,9=2,10=2,28=1/1;
* 1/52=2,53=1022/20;
* 3/6=3,11=9,16=1,25=1,30=1,116=-2/1,2,3;
* 4/5=1,17=6,18=6/1;
* 5/5=2,38=6/2;
* 6/7=2,8=2,9=2,10=2,28=1/1;
* 1/52=2,53=3014/20;
99/5=1,9=1,10=32/99;
</pre>

The asterisked lines can be removed as they correspond to the low real system and adjust the last line to read <nowiki>99/10=32/99</nowiki>. We can now construct the input file to punch out the high-model orbitals. '''Remember''' to add ''cp fort.7 $WORK/$FLD/$FLNM.orbs'' to your jobscript file after the gaussian execution line as this file will contain the punched orbitals.

<pre>
%nprocshared=2
%mem=2000MB
%chk=/work/lmt09/PHD_Y2/MALA_CYHEX/ONIOM/macyhexdiene_S0_SPpunch_oniom_cas66_sto3g_hf_sto3g
# nonstd
1/38=1,52=2/1;
2/12=2,15=1,17=6,18=5,40=1/2;
1/38=1,52=2,53=3172/20;
3/6=3,11=9,16=1,25=1,30=1,116=-2/1,2,3;
4/5=1,17=6,18=6/1;
5/5=2,38=6/2;
6/7=2,8=2,9=2,10=2,28=1/1;
1/52=2,53=2032/20;
3/6=3,16=1,25=1,32=1,116=101/1,2,3;
4/5=1,17=6,18=6/1,5;
5/5=2,17=1000000,38=6/10;
6/7=2,8=2,9=2,10=2,28=1/1;
99/10=32/99;

#p oniom(casscf(6,6)/sto-3g:hf/sto-3g) guess=read nosymm punch=MO
Punch high model orbitals for localization

0 1 0 1 0 1
H 0 -0.26330500 -1.99941700 -1.21363800 H
C 0 0.78714900 1.41335300 0.46807900 H
C 0 1.26328700 1.66455000 -0.88162300 H
O 0 -1.92502500 -0.01205900 2.50907200 L
C 0 2.10729300 0.78800800 -1.50119200 H
H 0 0.93300900 2.55923500 -1.39147400 H
C 0 2.36471400 -0.53645200 0.68140600 L H 9 0.0000
C 0 -0.44009500 -1.08924500 0.83298900 H
C 0 1.08572800 0.18779800 1.13172100 H
C 0 -0.78359900 -1.36617800 -0.51617700 H
C 0 2.60670600 -0.49404400 -0.84831600 L H 5 0.0000
O 0 -2.42068300 0.13236600 0.27717100 L
C 0 -1.63821400 -0.28747100 1.36303200 L H 8 0.0000
H 0 2.12198000 -1.34178900 -1.33290700 L
O 0 -2.47630400 -0.31974200 -1.96661200 L
H 0 3.67299000 -0.60346100 -1.04369100 L
C 0 -1.93186300 -0.50956600 -0.89422500 L H 10 0.0000
H 0 0.08918800 2.11204200 0.91605600 H
H 0 -0.01783200 -1.85015200 1.47955400 H
H 0 0.96794500 0.21509600 2.21057000 H
H 0 2.35971100 -1.57164700 1.01714300 L
H 0 3.19920800 -0.04865300 1.18673800 L
H 0 2.42869300 0.97055600 -2.52117800 H

</pre>

====Obtain Geometry====
Now we have the orbitals for the high model in the .orbs file, however, if we wish to visualize these with gaussview we need to have them in a Gaussian output and so we need the geometry of the high model system. This can be done by taking the geometry output by the optimization and using the ''onlyinputfiles'' option.
<pre>
%nprocshared=1
%mem=800MB
%chk=/work/lmt09/PHD_Y2/MALA_CYHEX/ONIOM/macyhexdiene_S0_SPinput_oniom_cas66_sto3g_hf_sto3g
#p oniom(casscf(6,6)/sto-3g:hf/sto-3g)=onlyinputfiles nosymm

Input files

0 1 0 1 0 1
H 0 -0.26330500 -1.99941700 -1.21363800 H
C 0 0.78714900 1.41335300 0.46807900 H
C 0 1.26328700 1.66455000 -0.88162300 H
O 0 -1.92502500 -0.01205900 2.50907200 L
C 0 2.10729300 0.78800800 -1.50119200 H
H 0 0.93300900 2.55923500 -1.39147400 H
C 0 2.36471400 -0.53645200 0.68140600 L H 9 0.0000
C 0 -0.44009500 -1.08924500 0.83298900 H
C 0 1.08572800 0.18779800 1.13172100 H
C 0 -0.78359900 -1.36617800 -0.51617700 H
C 0 2.60670600 -0.49404400 -0.84831600 L H 5 0.0000
O 0 -2.42068300 0.13236600 0.27717100 L
C 0 -1.63821400 -0.28747100 1.36303200 L H 8 0.0000
H 0 2.12198000 -1.34178900 -1.33290700 L
O 0 -2.47630400 -0.31974200 -1.96661200 L
H 0 3.67299000 -0.60346100 -1.04369100 L
C 0 -1.93186300 -0.50956600 -0.89422500 L H 10 0.0000
H 0 0.08918800 2.11204200 0.91605600 H
H 0 -0.01783200 -1.85015200 1.47955400 H
H 0 0.96794500 0.21509600 2.21057000 H
H 0 2.35971100 -1.57164700 1.01714300 L
H 0 3.19920800 -0.04865300 1.18673800 L
H 0 2.42869300 0.97055600 -2.52117800 H

</pre>

====Localizing Orbitals====
The final part is to combine these two calculations to produce a Gaussian output with the orbitals on the high model. The first step is to copy and paste the high model input file from the ''onlyinputfiles'' output. '''Remember''' this is a single layer calculation so ensure that no ONIOM keywords are present. Once this is done add ''guess=cards'' to the route (note that in the example below the IOps have been removed as they are not necessary) and copy the .orbs file to the bottom of the input file. If localized orbitals are desired the vaarious keywords can be added here. This results in the following input:
[[Media:Male_loc.gjf]] 

We can now visualize the results which reveal a ''p-orbital'' on each carbon of the high model region, showing the correct active space has bee chosen.
[[Media:Male_loc.log]]

File:Lepsgui.m

2014-02-11T01:24:13Z

Lmt09: uploaded a new version of "File:Lepsgui.m"

File:Lepsgui.fig

2014-02-11T01:23:26Z

Lmt09: uploaded a new version of "File:Lepsgui.fig"

Gaussian Bugs and Issues

2013-09-10T17:31:18Z

Lmt09: /* CASSCF */

==Gaussian Bugs and Issues==

Section to place any suspected bugs and errors

===CASSCF===

'''1.'''Quartic linear search updated in l103 seems to change which inital CASSCF guess is read (Lee).

Use IOp(4/5=1) to read from the checkpoint file. 
Only a problem with GDV versions (GDVH13 and GDVH23 checked), not G09 (G09A02 and G09C01 checked)

'''2.'''Freq and StateAverage kewords are incompatible. Fix needs to use Freq and IOp for SA (David).

'''3.'''NBO analysis with CAS - attempting to save NBOs to checkpoint file causes core dump.

Use IOp(3/32=5). 
Tested for gdv versions h01-h23 and g09 c01-d01.

===CIS===

'''1.'''Geomtry optimization gives a core dump when reading initial guess vectors from the rwf. (Lee)

Use IOp(9/49=3) to force read from the checkpoint file. 
The checkpoint file must be available to read from in the first step 
If this is not the case run a single point and then run the optimization using the SP chk file with the IOp 
Checked with GDVH13 only.

===Opt===

====Modredundant====

'''1.''' Using the Frozen coordinate command whilst setting the value for the coordinate (i.e. B 1 2 2.2 F which should freeze the bond distance between atoms 1 and 2 to 2.2 angstrom) sets the correct angle but does not freeze it. (Lee)

Build the coordinate and freeze in different specifications e.g. 
B 1 2 2.2 B 
B 1 2 F

Computing Transition States in ONIOM

2013-02-20T18:21:47Z

Lmt09:

This documents list the method I use for finding TSs with ONIOM.

1. Run a scan over the TS coordinate.<br\>
If you get any discontinuities instead of a smooth curve around where the TS should be then you are scanning along the wrong coordinate
Look for atoms that move a large amount in the discontinuous step and scan along this one instead.

2. Take the step that has the highest energy and the lowest RMS force (the most TS like structure) and optimize this with the TS coordinate frozen.
Use geom=modredundant and add the specification under the molecule geometries.

3. Numerical compute the Hessian around the TS coordinate using geom=modredundant and specifying it under the coordinates.
Use Opt=(ts,noeigen), also things to consider are the step size and the use of quadmac which seems to be pretty much required

So basically it is the same as finding a normal TS but your options are a bit more restricted (e.g. no QST2).

Computing Transition States in ONIOM

2013-02-20T18:20:30Z

Lmt09:

Computing Transition States in ONIOM

2013-02-20T18:17:26Z

Lmt09: Created page with "This documents list the method I use for finding TSs with ONIOM. 1. Run a scan over the TS coordinate.<br\> If you get any discontinuities instead of a smooth curve around where..."

ONIOM

2013-02-20T18:12:00Z

Lmt09: /* Tutorials */

==Background Information==
[[Introduction]]<br\>
[[ONIOM for excited states]]<br\>
[[ONIOM for crossings]]<br\>

==Tutorials==
[https://www.ch.ic.ac.uk/wiki/index.php/First_steps_with_ONIOM First Steps With ONIOM]<br\>
[https://www.ch.ic.ac.uk/wiki/index.php/First_steps_with_ONIOM_:_excited_state_of_Bicyclo%2810%2C2%2C2%29Hexadeca-1%2815%29%2C12%2816%2913-triene Study of bicyclo(10,2,2)hexadeca-1(15),12(16)13-triene]<br\>
[https://wiki.ch.ic.ac.uk/wiki/index.php?title=Let%27s_try_with_a_smaller_cycle Study of bicyclo(8,2,2)hexadeca-1(13),10(14)11-triene]<br\>
[[Visualizing High Model Orbitals]]<br\>
[[Computing Transition States in ONIOM]]

== ONIOM for Biomolecules ==

===Creating Input Files===
[[Guide to Creating ONIOM input files for biomolecules]]

===Force Field Methods===
[[AMBER]]

===Systems (older notes)===
[[PYP : Photoactive yellow protein]] 
[[GFP : Green Fluorescent Protein]]

==ONIOM Frequency Calculations==
[[Frequency Calculations on a Small System (high level vs. ONIOM)]]

==Know Problems==
[[gdvh01_pgi_725]] 
[[gdvh08_pgi_806]]

Mod:phys3

2013-01-09T15:51:13Z

Lmt09: /* Optimizing the "Chair" and "Boat" Transition Structures */

See also:[[Mod:timetable|Timetable]], [[Mod:lectures|Intro lecture]], [[mod:programs|Programs]], [[mod:organic|Module 1]], [[Mod:inorganic|Module 2]], [[Mod:phys3|Module 3]], [http://www.gaussian.com/g_tech/gv5ref/gv5ref_toc.htm Gaussian Online User Manual] | [http://faculty.ycp.edu/~jforesma/educ/visual/index.html Visualization Tutorials]
= Module 3 =

In this set of computational experiments, you will characterise transition structures on potential energy surfaces for the Cope rearrangement and Diels Alder cycloaddition reactions.

There are two parts:
a) tutorial material: how to use the programs and methods,
b) more challenging examples, with guidelines but fewer explicit instructions.



In the second year physical chemistry laboratory, you may have carried out dynamics calculations using model potential energy surfaces to explore transition states. In that computational experiment, the total energy could quickly be calculated for different geometries of a triatomic system using an analytical function of the atomic coordinates (for more information, see for example [http://books.google.com/books?id=T8IZ1aa_FRkC&pg=RA1-PA36&lpg=RA1-PA36&dq=%22lake+eyring%22&source=web&ots=OXY00lSZ7D&sig=Ld_MTNwNjUDNGzB_5w1IxaMBMPU&hl=en&sa=X&oi=book_result&resnum=7&ct=result here] and [http://www.rsc.org/ejarchive/DC/1979/DC9796700007.pdf here]).

In this experiment, you will be studying transition structures in larger molecules. There are no longer fitted formulae for the energy, and the molecular mechanics / force field methods that work well for structure determination cannot be used (in general) as they do not describe bonds being made and broken, and changes in bonding type / electron distribution. (This is the main difference from Module 1). Instead, we use molecular orbital-based methods, numerically solving the Schrodinger equation, and locating transition structures based on the local shape of a potential energy surface. As well as showing what transition structures look like, reaction paths and barrier heights can also be calculated.

==The Cope Rearrangement Tutorial==



''This part of the module is described as a 'tutorial' because it's an introduction to various computational techniques for locating transition structures on potential energy surfaces. It's different to the GaussView tutorial you may have worked through earlier: it's an exercise where you're given specific instructions, see if you can follow them, and also whether there are problems or better ways of carrying the exercise out. Please include this part in your write-up. Marks will be given for correct answers, the documentation showing how you got these, discussion, and how you went about solving any problems you encountered.''

In this tutorial we will use the Cope rearrangement of 1,5-hexadiene as an example of how to study a chemical reactivity problem.

Your objectives are to locate the low-energy minima and transition structures on the C6H10 potential energy surface, to determine the preferred reaction mechanism.

[[Image:pic1.jpg|right|thumb|Cope rearrangement]]

This [3,3]-sigmatropic shift rearrangement has been the subject of numerous experimental and computational studies (e.g. Houk et al. {{DOI|10.1021/ja00101a078}}), and for a long time its mechanism (concerted, stepwise or dissociative) was the subject of some controversy. Nowadays it is generally accepted that the reaction occurs in a concerted fashion via either a "chair" or a "boat" transition structure, with the "boat" transition structure lying several kcal/mol higher in energy. The B3LYP/6-31G* level of theory has been shown to give activation energies and enthalpies in remarkably good agreement with experiment. In this tutorial we will show how these can be calculated using Gaussian.

 

{| align="center"
|- align="center"
|
[[Image:pic2a.jpg]]
|
[[Image:pic2b.jpg]]
|- align="center"
| align="center" | ''Chair Transition State''
| align="center" | ''Boat Transition State''
|}

 

===Optimizing the Reactants and Products===

'' In this section you will learn how to optimize a structure, symmetrize it to find its point group, calculate and visualize vibrational frequencies and correct potential energies in order to compare them with experimental values. It is assumed that you are already familiar with using the builder in GaussView. ''

 

(a) Using GaussView, draw a molecule of 1,5-hexadiene with an "anti" linkage (aproximately a.p.p conformation) for the central four C atoms . Clean the structure using the '''Clean''' function under the '''Edit''' menu.

Now we will optimize the structure at the '''HF/3-21G''' level of theory. Select '''Gaussian''' under the '''Calculate''' menu, click on the '''Job Type''' tab and choose '''Optimization'''. The default method should already be Hartree Fock and the default basis set is 3-21G, so there should be no need to change these. You can check this by clicking on the '''Method''' tab. Change the %mem under the '''Link 0''' tab to 250 MB (though this could be increased to 500 MB). Submit the job by clicking on the '''Submit''' button at the bottom of the window and give the job a meaningful name (e.g. react_anti).

When the job has finished, you will be asked if you want to open a file. Select '''Yes''' and choose the checkpoint (chk) file with the name of the job you have just run (e.g. react_anti.chk). This checkpoint file is a binary file that stores data calculated by Gaussian. The name of the chk file should have been assigned by default, but by default, this file will be created in the C:\Windows\G03\Scratch folder. Once the file has been opened, click on the '''Summary''' button under the '''Results''' menu and make a note of the energy.


Does your final structure have symmetry? Select '''Symmetrize''' under the '''Edit''' menu (note that sometimes it is necessary to relax the search criteria under the '''Point Group''' menu). Make a note of the point group.

 

(b) Now draw another molecule of 1,5-hexadiene with a "gauche" linkage for the central four C atoms. Would you expect this structure to have a lower or a higher energy than the anti structure you have just optimized? Optimize the structure at the '''HF/3-21G''' level of theory and compare your final energy with that obtained in (a). Again, check if the molecule has symmetry and make a note of the point group.

 

(c) Normally, calculated activation energies and enthalpies use the lowest energy conformation of a reactant molecule as a reference. Based on your results from above, try to predict what the lowest energy conformation of 1,5-hexadiene might be. Test out your hypothesis by drawing the structure and optimizing it.

 

(d) A table containing the low energy conformers of 1,5-hexadiene and their point groups is shown in [[Mod:phys3#Appendix 1|Appendix 1]]. Compare the structures that you have optimized with those in the table and see if you can identify your structure.

 

(e) Draw the Ci ''anti2'' conformation of 1,5-hexadiene (unless you have already located it). Optimize it at the '''HF/3-21G''' level of theory and make sure it has Ci symmetry. Compare your final energy to the one given in the table.


 

(f) When you are happy that your structure is the same as the one in the table, reoptimize it at the '''B3LYP/6-31G*''' level (6-31G* is equivalent to 6-31G(d) by selecting '''DFT''' under the '''Method''' menu and '''B3LYP''' from the box with the functionals on the right-hand side. Now select '''Link 0''' and change the name of the chk file to the name of the DFT optimization that you are about to run. Note that it is always advisable to do this when re-using or modifying existing structures to ensure that the original chk file is not overwritten. Run the job and make a note of the energy. Now compare the final structures from the '''HF/3-21G''' calculation with that at the higher level of theory. How much does the overall geometry change?

 

(g) The final energies given in the output file represent the energy of the molecule on the bare potential energy surface. To be able to compare these energies with experimentally measured quantities, they need to include some additional terms, which requires a frequency calculation to be carried out. The frequency calculation can also be used to characterize the critical point, i.e. to confirm that it is a minimum in this case: that all vibrational frequencies are real and positive.

Starting from your optimized B3LYP/6-31G* structure, run a frequency calculation at the same level of theory. You can do this by selecting '''Frequency''' under the '''Job Type''' tab. Ensure that the method is still correctly specified under the '''Method''' tab (''caution: on Windows, sometimes 'scrf=(solvent=water,check)' is incorrectly added!'') and then change the name of the chk file under the '''Link 0''' tab to the name of the frequency job that you are about to run. Run the job. Once the job has finished, open the log file this time. Select '''Vibrations''' under the '''Results''' menu. A list of all the vibrational frequencies modes should appear. Check that there are no imaginary frequencies, only real ones. You can visualize some of these vibrations under this menu and simulate the infrared spectrum.


Now, select '''View File''' under the '''Results''' menu and open the output file in the visualizer. Scroll down to the section beginning '''Thermochemistry'''. Under the vibrational temperatures a list of energies should be printed. Make a note of (i) the sum of electronic and zero-point energies, (ii) the sum of electronic and thermal energies, (iii) the sum of electronic and thermal enthalpies, and (iv) the sum of electronic and thermal free energies. The first of these is the potential energy at 0 K including the zero-point vibrational energy (E = Eelec + ZPE), the second is the energy at 298.15 K and 1 atm of pressure which includes contributions from the translational, rotational, and vibrational energy modes at this temperature (E = E + Evib + Erot + Etrans), the third contains an additional correction for RT (H = E + RT) which is particularly important when looking at dissociation reactions, and the last includes the entropic contribution to the free energy (G = H - TS). It is important to make sure that you select the correct energy/enthalpy term to compare to your experimental values. Note that these corrections can also be calculated at other temperatures using the '''Temperature''' option in Gaussian, If you have time, try re-calculate these quantities at 0 K as shown in the [[mod:gv_advanced | Advanced GaussView Tutorial]].

 

===Optimizing the "Chair" and "Boat" Transition Structures ===

'' In this section you will learn how to set up a transition structure optimization (i) by computing the force constants at the beginning of the calculation, (ii) using the redundant coordinate editor, and (iii) using QST2. You will also visualize the reaction coordinate and run the IRC (Intrinisic Reaction Coordinate) and calculate the activation energies for the Cope rearrangement via the "chair" and "boat" transition structures. ''

 

The "chair" and "boat" transition structures for the Cope rearrangement are shown in [[Mod:phys3#Appendix 2|Appendix 2]]. Both consist of two C3H5 allyl fragments positioned approximately 2.2 Å apart, one with C2h symmetry and the other with C2v symmetry.

 

(a) Draw an allyl fragment (CH2CHCH2) and optimize it using the '''HF/3-21G''' level of theory. Your structure should look like one half of the transition structures shown below.

Now open a new GaussView window by going to the '''File''' menu and selecting '''New''' and then '''Create MolGroup'''. Copy the optimized allyl structure from the first calculation by selecting '''Copy''' under the '''Edit''' menu, and then paste it twice into the new window by selecting '''Paste''' and then '''Append Molecule'''. Now orient the two fragments so that they look roughly like the chair transition state below by using the '''Shift Alt keys + Left Mouse button''' to translate one fragment with respect to the other and the '''Alt key + Left Mouse button''' to rotate it. The distance between the terminal ends of the allyl fragments should be approximately 2.2 Å apart. Save this structure to a Gaussian input file with a meaningful name (e.g. chair_ts_guess).

We are now going to optimize this transition state manually in two different ways. Transition state optimizations are more difficult than minimizations because the calculation needs to know where the negative direction of curvature (i.e. the reaction coordinate) is. If you have a reasonable guess for your transition structure geometry, then normally the easiest way to produce this information is to compute the force constant matrix (also known as the Hessian) in the first step of the optimization which will then be updated as the optimization proceeds. This is what we will try to do in the next section. However, if the guess structure for the transition structure is far from the exact structure, then this approach may not work as the curvature of the surface may be significantly different at points far removed from the transition structure. In some cases, a better transition structure can be generated by freezing the reaction coordinate (using '''Opt=ModRedundant''' and minimizing the rest of the molecule. Once the molecule is fully relaxed, the reaction coordinate can then be unfrozen and the transition state optimization is started again. One advantage of doing this, is that it may not be necessary to compute the whole Hessian once this has been done, and just differentiating along the reaction coordinate might give a good enough guess for the initial force constant matrix. This can save a considerable amount of time in cases where the force constant calculation is expensive.

 

(b) Use Hartree Fock and the default basis set 3-21G for parts (b) to (f).

Create a new MolGroup ('''File → New → Create MolGroup''') and copy and paste your guess structure into the window. Now set up a Gaussian optimization for a transition state. Go to the '''Gaussian''' menu under '''Calculate''' and click on the '''Job Type''' tab. Select '''Opt+Freq''' and then change '''Optimization to a Minimum''' to '''Optimization to a TS (Berny)'''. Choose to calculate the force constants '''Once''' and in the Additional keyword box at the bottom, type '''Opt=NoEigen'''. The latter stops the calculation crashing if more than one imaginary frequency is detected during the optimization which can often happen if the guess transition structure is not good enough. Submit the job. If the job completes successfully, you should have optimized to the structure shown in [[Mod:phys3#Appendix 2|Appendix 2]] and the frequency calculation should give an imaginary frequency of magnitude 818 cm-1. Animate the vibration and ensure that it is the one corresponding to the Cope rearrangement.

 

(c) Now we will try optimizing the transition structure again using the frozen coordinate method. Create a new MolGroup ('''File → New → Create MolGroup''') and copy and paste your guess structure into the window again. Now select '''Redundant Coord Editor''' from the '''Edit''' menu. Click on the highlighted file icon at the top left-hand corner (Create a New Coordinate) and a line should appear below saying '''Add Unidentified (?, ?, ?, ?)'''. Now go back to the GaussView window and select two of the terminal carbons from the allyl fragments which form/break a bond during the rearrangement. Return to the coordinate editor and select '''Bond''' instead of '''Unidentified''' and select '''Freeze Coordinate''' instead of '''Add'''. Now click on the icon again to generate another coordinate. This time select the opposite two terminal atoms and again select '''Bond''' and '''Freeze Coordinate'''. Click OK. Now set up the optimization as if it were a minimum and you should see the option '''Opt=ModRedundant''' already included in the input line. Submit the job.

'''Note:''' GaussView allows you to produce an input file with the frozen coordinate specified as e.g. <tt>B 5 1 2.200000 F</tt>. Unfortunately, a recent update to the Gaussian program means it does not recognise this syntax, and just ignores this line. This means that the coordinate ends up being optimised rather than frozen. Therefore do not use this method, but ensure the guess structure has suitable guess transition bond distances(~2.2 Å) using the ''Modify Bond'' tool in GaussView --[[User:Rzepa|Rzepa]] 14:39, 29 October 2012 (UTC)

'''Note 2:''' If you set the coordinate distance (i.e. B 5 1 2.2 B) and freeze it (i.e. B 5 1 F) as separate inputs to the modredundant editor, it appears to work as expected. (Lee)

 

(d) When the job has finished, open the chk file. You should find that the optimized structure looks a lot like the transition you optimized in section (b), except the bond forming/breaking distances are fixed to 2.2 Å. Now we are going to optimize them too. Open the '''Redundant Coord Editor''' from the '''Edit''' menu again and create a new coordinate as before by clicking on the icon, Select one of the bonds that was previously frozen and this time choose '''Bond''' instead of '''Unidentified''' and '''Derivative''' instead of '''Add'''. Repeat the procedure for the other bond. This time you need to set up a transition state optimization but we are not going to calculate the force constants as we did in section (b) (so we leave this option as '''Never'''), instead we will use a normal guess Hessian modified to include the information about the two coordinates we are differentiating along. Change the name of the chk file in '''Link 0''' if you do not want to write over the previous calculation and submit the job. When the calculation has finished, open the chk file, check the bond forming/bond breaking bond lengths and compare the structure to the one you optimized in section (b).

 

(e) Now we will optimize the boat transition structure. We will do this using the '''QST2''' method. In this method, you can specify the reactants and products for a reaction and the calculation will interpolate between the two structures to try to find the transition state between them. You must make sure that your reactants and products are numbered in the same way. Therefore, although our reactants and products are both 1,5-hexadiene, we will need to manually change the numbering for the product molecule so that it corresponds to the numbering obtained if our reactant had rearranged.

e.g.

<center>[[Image:pic3.jpg|200px]]</center>

Open the chk file corresponding to the optimized Ci reactant molecule (''anti2'' in [[Mod:phys3#Appendix 1|Appendix 1]]). Now open a second window and create a new MolGroup. Copy the optimized reactant molecule into the new window. In the same window, now select '''File → New → Add to MolGroup'''. The original molecule should disappear and a green circle should appear at the top left-hand corner with a '''2''' next to it. Clicking on the down arrow by the '''2''' will take you back to the original window and you will see your molecule again. This is how we read multiple geometries into GaussView. Go back to window '''2''', and copy and paste the reactant molecule a second time. This is going to be the product molecule and will be the molecule on which we need to change the numbering. If you now click on the icon showing two molecules side by side, then you can view both molecules simultaneously.

Now go to the '''View''' menu and select '''Labels''' so that you can see the numbering on both structures. Orient the two structures separately so they look something like the following:

{| align="center"
|- align="center"
|
[[Image:pic4a.jpg|200px]]
|
[[Image:pic4b.jpg|200px]]
|- align="center"
| align="center" | ''Reactant''
| align="center" | ''Product''
|}

 

Now click on the product structure. Go to the '''Edit''' menu and select '''Atom List'''. Starting from Atom 1 on the reactant, go through and renumber all the atoms on the Product so that they match the reactant molecule, e.g. for the numbering above you would start by changing atom '''6''' on the product molecule to atom '''3'''. The other atom numbers will update as you do this so make sure you do it in the correct order. At the end, the numbering on your two molecules should correspond to each other in the following way:

 

{| align="center"
|- align="center"
|
[[Image:pic5a.jpg|200px]]
|
[[Image:pic5b.jpg|200px]]
|- align="center"
| align="center" | ''Reactant''
| align="center" | ''Product''
|}

 

Now we will set up the first '''QST2''' calculation. Go to the '''Gaussian''' menu and select '''Job Type''' as '''Opt+Freq''', and optimize to a transition state. This time you will have two options - '''TS (Berny)''' which we used in the previous calculations and '''TS (QST2)'''. Select '''TS (QST2)'''. Submit the job.

You will find that the job fails. To see why, open the chk file you created and view the structure. You will see that it looks a bit like the chair transition structure but more dissociated. In fact when the calculation linearly interpolated between the two structures, it simply translated the top '''allyl''' fragment and did not even consider the possibility of a rotation around the central bonds. It is clear that the QST2 method is never going to locate the boat transition structure if we start from these reactant and product structures.

Now go back to the original input file where you set up your QST2 calculation. We will now modify the reactant and product geometries so that they are closer to the boat transition structure. Click on the reactant molecule first and select the central '''C-C-C-C''' dihedral angle (i.e. '''C2-C3-C4-C5''' for the molecule above) and change the angle to 0o. Then select the inside '''C-C-C''' (i.e. '''C2-C3-C4''' and '''C3-C4-C5''' for the molecule above) and reduce them to 100o. Do the same for the product molecule. Your reactant and product molecules should now look like the following:

 

{| align="center"
|- align="center"
|
[[Image:pic6a.jpg|200px]]
|
[[Image:pic6b.jpg|200px]]
|- align="center"
| align="center" | ''Reactant''
| align="center" | ''Product''
|}

 

Set up the QST2 calculation again, renaming both the chk file under '''Link 0''' and the input file. Run the job again. This time it should converge to the boat transition structure. Check that there is only one imaginary frequency and visualize its motion.

The object of this exercise is to illustrate that although the QST2 method is has some advantages because it is fully automated, it can often fail if your reactants and products are not close to the transition structure. There is another method, the '''QST3''' method, that allows you to input the geometry of a guess transition structure also and this can often be more reliable. If you have time, you can try generating a guess boat transition structure and see if you can get the calculation to converge using the original reactant and product molecules. Remember to check the atom numbers in the transition structure are in the right order.

 

(f) Take a look at your optimized chair and boat transition structures. Which conformers of 1,5-hexadiene do you think they connect? You will find that it is almost impossible to predict which conformer the reaction paths from the transitions structures will lead to. However, there is a method implemented in Gaussian which allows you to follow the minimum energy path from a transition structure down to its local minimum on a potential energy surface. This is called the '''Intrinsic Reaction Coordinate''' or '''IRC''' method. This creates a series of points by taking small geometry steps in the direction where the gradient or slope of the energy surface is steepest.

Open the chk file for one of your optimized chair transition structures. Under the '''Gaussian''' menu, select '''IRC''' under the '''Job Type''' tab. You will be presented with a number of options. The first is to decide whether to compute the reaction coordinate in one or both directions. As our reaction coordinate is symmetrical, we will only choose to compute it in the forward direction. Normally you would do both forward and reverse, either in one job or in two separate jobs. You are also given the option to calculate the force constants once, at every step along the IRC or to read them from the chk file, in this case choose calculate always. You would use the latter option if you have previously run a frequency calculation. (The '''IRCMax''' option can also be specified here. This takes a transition structure as its input, and finds the maximum energy along a specified reaction path, taking into account zero-point energy etc., and produces all the quantities needed for a variational transition state theory calculation. We will leave this unchecked for the purposes of this exercise.) The final option to consider is the number of points along the IRC. The default is '''6''' but this is normally never enough. Let's change this to 50 and see how the calculation progresses. Change the name of the chk file under '''Link 0''' and submit the job. The job will take a while so now is a good time to take a coffee break...

When the IRC calculation has finished, open the chk file with all the intermediate geometries and see how the calculation has progressed. You will find that it hasn't reached a minimum geometry yet. This leaves you three options: (i) you can take the last point on the IRC and run a normal minimization; (ii) you can restart the IRC and specify a larger number of points until it reaches a minimum; (iii) you can redo the IRC specifying that you want to compute the force constants at every step. There are advantages and disadvantages to each of these approaches. Approach (i) is the fastest, but if you are not close enough to a local minimum, you may end up in the wrong minimum. Approach (ii) is more reliable but if too many points are needed, then you can also veer off in the wrong direction after a while and end up at the wrong structure. Approach (iii) is the most reliable but also the most expensive and is not always feasible for large systems. You can try any or all of these approaches and see which conformation you end up in.

 

(g) Finally we need to calculate the activation energies for our reaction via both transition structures. To do this we will need to reoptimize the chair and boat transition structures using the '''B3LYP/6-31G*''' level of theory and to carry out frequency calculations. You can start from the HF/3-21G optimized structures. Once the calculations have converged, compare both the geometries and the difference in energies between the reactants and transition states at the two levels of theory. What you should find is that the geometries are reasonably similar, but the energy differences are markedly different.

As a consequence of this, it is often more computational efficient to map the potential energy surface using the low level of theory first and then to reoptimize at the higher level as we have done in this exercise.

The experimental activation energies are 33.5 ± 0.5 kcal/mol via the chair transition structure and 44.7 ± 2.0 kcal/mol via the boat transition structure at 0 K. If you take the values computed at 0 K, how close are they to the experimental values? You can also find the energies with thermal correction at 298.15 K under the Thermochemistry data in the output file. If you have time, you can recompute them at higher temperature. Alternatively, you can use the utility program '''FreqChk''' to obtain energies at a different temperature. This only requires the chk file from a frequency calculation and allows you to retrieve frequency and thermochemistry data as well as calculating them with an alternate temperature, pressure, scale factor, and/or isotope substitutions. The '''FreqChk''' utility program can be accessed from '''Gaussian03W'''. Launch '''Gaussian03W'''. Select '''utilities''' from the menu and click on '''FreqChk''' to launch the utility program. You will be prompted for a chk file. Select your chk file from the C:\G03W\Scratch directory and follow the instructions from this [http://www.gaussian.com/g_tech/g_ur/u_freqchk.htm web link] to proceed.

 

=== Appendix 1 ===

 

{| border="1" cellpadding="5"
|- align="center"
| width="150" | '''Conformer'''
| width="150" | '''Structure'''
| width="100" | '''Point Group'''
| width="200" | '''Energy/Hartrees HF/3-21G'''
| width="200" | '''Relative Energy/kcal/mol'''
|- align="center"
| ''gauche1''
|
[[Image:gauche1.jpg|150px]]
| C2
| -231.68772
| 3.10
|- align="center"
| ''gauche2''
|
[[Image:gauche2.jpg|150px]]
| C2
| -231.69167
| 0.62
|- align="center"
| ''gauche3''
|
[[Image:gauche3.jpg|150px]]
| C1
| -231.69266
| 0.00
|- align="center"
| ''gauche4''
|
[[Image:gauche4.jpg|150px]]
| C2
| -231.69153
| 0.71
|- align="center"
| ''gauche5''
|
[[Image:gauche5.jpg|150px]]
| C1
| -231.68962
| 1.91
|- align="center"
| ''gauche6''
|
[[Image:gauche6.jpg|150px]]
| C1
| -231.68916
| 2.20
|- align="center"
| ''anti1''
|
[[Image:anti1.jpg|150px]]
| C2
| -231.69260
| 0.04
|- align="center"
| ''anti2''
|
[[Image:anti2.jpg|150px]]
| Ci
| -231.69254
| 0.08
|- align="center"
| ''anti3''
|
[[Image:anti3.jpg|150px]]
| C2h
| -231.68907
| 2.25
|- align="center"
| ''anti4''
|
[[Image:anti4.jpg|150px]]
| C1
| -231.69097
| 1.06
|}

 

=== Appendix 2 ===

{| cellpadding="5"
|- align="center"
|
[[Image:appendix2a.jpg|300px]]
|- align="center"
| ''C2h Chair Transition State''
|- align="center"
|- align="center"
|- align="center"
|
[[Image:appendix2b.jpg|300px]]
|- align="center"
| ''C2v Boat Transition State''
|}

 

=== Results Table ===

 

''' Summary of energies (in hartree) '''

{| border="1" cellspacing="1" cellpadding="10"
|-
| ''' '''
!colspan="3"|'''HF/3-21G'''
!colspan="3"|'''B3LYP/6-31G*'''
|-
| ''' '''
| width="125" align="center" | '''Electronic energy'''
| width="125" align="center" | '''Sum of electronic and zero-point energies'''
| width="125" align="center" | '''Sum of electronic and thermal energies'''
| width="125" align="center" | '''Electronic energy'''
| width="125" align="center" | '''Sum of electronic and zero-point energies'''
| width="125" align="center" | '''Sum of electronic and thermal energies'''
|-
| ''' '''
| ''' '''
| width="125" align="center" | '''at 0 K'''
| width="125" align="center" | '''at 298.15 K'''
| ''' '''
| width="125" align="center" | '''at 0 K'''
| width="125" align="center" | '''at 298.15 K'''
|-
| '''Chair TS'''
| align="center" | -231.619322
| align="center" | -231.466705
| align="center" | -231.461346
| align="center" | -234.556983
| align="center" | -234.414919
| align="center" | -234.408998
|-
| '''Boat TS'''
| align="center" | -231.602802
| align="center" | -231.450929
| align="center" | -231.445300
| align="center" | -234.543093
| align="center" | -234.402340
| align="center" | -234.396006
|-
| '''Reactant (''anti2'')'''
| align="center" | -231.692535
| align="center" | -231.539539
| align="center" | -231.532566
| align="center" | -234.611710
| align="center" | -234.469203
| align="center" | -234.461856
|}

 *1 hartree = 627.509 kcal/mol 

''' Summary of activation energies (in kcal/mol) '''

{| border="1" cellspacing="1" cellpadding="10"
|-
| ''' '''
| width="125" align="center" | '''HF/3-21G'''
| width="125" align="center" | '''HF/3-21G'''
| width="125" align="center" | '''B3LYP/6-31G*'''
| width="125" align="center" | '''B3LYP/6-31G*'''
| width="125" align="center" | '''Expt.'''
|-
| ''' '''
| width="125" align="center" | '''at 0 K '''
| width="125" align="center" | '''at 298.15 K'''
| width="125" align="center" | '''at 0 K'''
| width="125" align="center" | '''at 298.15 K'''
| width="125" align="center" | '''at 0 K'''
|-
| '''ΔE (Chair)'''
| align="center" | 45.70
| align="center" | 44.69
| align="center" | 34.06
| align="center" | 33.17
| align="center" | 33.5 ± 0.5
|-
| '''ΔE (Boat)'''
| align="center" | 55.60
| align="center" | 54.76
| align="center" | 41.96
| align="center" | 41.32
| align="center" | 44.7 ± 2.0
|}

 

==The Diels Alder Cycloaddition==

In this exercise, you will characterise transition structures using any of the methods described above in the tutorial: the choice is up to you. In addition, you will look at the shape of some of the molecular orbitals. To help you structure your report, there is a data/discussion sheet at the end of this section.

[[Image:mb_da1.jpg |right|thumb|Diels Alder cycloaddition]]
The Diels Alder reaction belongs to a class of reactions known as pericyclic reactions. The π orbitals of the dieneophile are used to form new σ bonds with the π orbitals of the diene. Whether or not the reactions occur in a concerted stereospecific fashion ('''allowed''') or not ('''forbidden''') depends on the number of π electrons involved. In general the HOMO/LUMO of one fragment interacts with the HOMO/LUMO of the other reactant to form two new bonding and anti-bonding MOs. The nodal properties allow one to make predictions according to the following rule:

''If the HOMO of one reactant can interact with the LUMO of the other reactant then the reaction is '''allowed'''.''

''The HOMO-LUMO can only interact when there is a significant overlap density. If the orbitals have different symmetry properties then no overlap density is possible and the reaction is '''forbidden'''.''

If the dieneophile is substituted, with substituents that have π orbitals that can interact with the new double bond that is being formed in the product, then this interaction can stabilise the regiochemistry (i.e. head to tail versus tail to head) of the reaction. In this exercise you will study the nature of the transition structure for the Diels Alder reaction, both for the prototypical reaction and for the case where both diene and dieneophile carry substituents, and where secondary orbital effects are possible. Clearly, the factors that control the nature of the transition state are quantum mechanical in origin and thus we shall use methods based upon quantum chemistry.

 

Shown on the right is a diagram of the transition state for the Diels-Alder reaction between ethylene and butadiene. The ethylene approaches the cis form of butadiene from above.
[[Image:mb_da2.jpg |right|thumb|Ethylene+Butadiene cycloaddition]]

Before beginning our quantitative study, it is helpful to discuss the interaction of the π orbitals in a simple qualitative way.

''You will confirm some of these considerations in your computations.''

The principal orbital interactions involve the π/ π* orbitals of ethylene and the HOMO/LUMO of butadiene. It is referred to as [4s + 2s] since one has 4 π orbitals in the π system of butadiene. The orbitals of ethylene and butadiene and ethylene can be classified as symmetric '''s''' or anti-symmetric '''a''' with respect to the plane of symmetry shown.

The HOMO of ethylene and the LUMO of butadiene are both '''s''' (symmetric with respect to the reflection plane) and the LUMO of ethylene and the HOMO of butadiene are both '''a'''. Thus it is the HOMO-LUMO pairs of orbital that interact, and energetically, the HOMO of the resulting adduct with two new σ bonds is '''a'''.

 

=== Exercise ===

Use the the AM1 semi-empirical molecular orbital method for these calculations (to start with).

i) Use GaussView to build cis butadiene, and optimize the geometry using Gaussian. Plot the HOMO and LUMO of cis butadiene and determine its symmetry (symmetric or anti-symmetric) with respect to plane.

''There are two ways to do this in GaussView. One is: Select '''Edit→MOs'''. Select the HOMO and the LUMO from the MO list (highlights it yellow). Click the button '''Visualise''' (not Calculation), then '''Update'''. Alternately, having calculated the surface for this orbital, you can display it in the main GaussView window for the molecule, from the '''Results→Surfaces''' menu. Select '''Surface Actions→Show Surface'''. Having displayed the surface this way, you can also select '''View→Display Format→Surface''', and change '''Solid''' to '''Mesh'''.''

 

ii) Computation of the Transition State geometry for the prototype reaction and an examination of the nature of the reaction path.

[[Image:mb_da3.jpg |right|thumb|]]

The transition structure has an envelope type structure, which maximizes the overlap between the ethylene π orbitals and the π system of butadiene. One way to obtain the starting geometry is to build the bicyclo system (b) and then remove the -CH2-CH2- fragment. One must then guess the interfragment distance (dashed lines) and optimize the structure, but use any method you wish, based on the tutorial above, to characterise the transition structure. Confirm you have obtained a transition structure for the Diels Alder reaction!

[[Image:mb_da4.jpg |right|thumb|guessing the transition structure]]

Once you have obtained the correct structure, plot the HOMO as in (i). Rotate the molecule so that the symmetry and nodal properties of the system can be interpreted, and save a copy of the image.

 

(iii) To Study the regioselectivity of the Diels Alder Reaction

Cyclohexa-1,3-diene '''1''' undergoes facile reaction with maleic anhydride '''2''' to give primarily the endo adduct. The reaction is supposed to be kinetically controlled so that the exo transition state should be higher in energy.

[[Image:Bearpark_pic_edit_by_jm906.JPG |right|thumb|regioslectivity]]

Locate the transition structures for both 3 and 4. Compare the energies of the endo and exo forms.

Measure the bond lengths of the partly formed σ C-C bonds and the other C-C distances. Make a sketch with the important bond lengths. Measure the orientation, (C-C through space distances between the -(C=O)-O-(C=O)- fragment of the maleic anhydride and the C atoms of the “opposite” -CH2-CH2- for the exo and the “opposite” -CH=CH- for the endo). The structure must be a compromise between steric repulsions of the -CH2-CH2- fragment and the maleic anhydride for the exo versus secondary orbital interactions between the π systems of -CH=CH- and -(C=O)-O-(C=O)- fragment for the endo.

Plot the HOMO as in the previous exercise. Examine carefully the nodal properties of the HOMO between the -(C=O)-O-(C=O)- fragment and the remainder of the system. What can you conclude about the so called “secondary orbital overlap effect”?

=== Suggested Discussion ===

Use this template as a guide. Screen images can be saved from the GaussView '''File''' menu.

''For cis butadiene'': 
Plot the HOMO and LUMO and determine the symmetry (symmetric or anti-symmetric) with respect to the plane.

''For the ethylene+cis butadiene transition structure'': 
Sketch HOMO and LUMO, labeling each as symmetric or anti symmetric.

Show the geometry of the transition structure, including the bond-lengths of the partly formed σ C-C bonds.

What are typical sp3 and sp2 C-C bondlengths? What is the van der Waals radius of the C atom? What can you conclude about the C-C bondlength of the partly formed σ C-C bonds in the TS.

Illustrate the vibration that corresponds to the reaction path at the transition state.
Is the formation of the two bonds synchronous or asynchronous?
How does this compare with the lowest positive frequency?

Is the HOMO at the transition structure '''s''' or '''a'''?

Which MOs of butadiene and ethylene have been used to form this MO?
Explain why the reaction is allowed.

''For the cyclohexa-1,3-diene reaction with maleic anhydride'': 
Give the relative energies of the exo and endo transition structures.
Comment on the structural difference between the endo and exo form. Why do you think that the exo form could be more strained?
Examine carefully the nodal properties of the HOMO between the -(C=O)-O-(C=O)- fragment and the remainder of the system. What can you conclude about the so called “secondary orbital overlap effect”?
(There is some discussion of this in Ian Fleming's book 'Frontier Orbitals and Organic Chemical Reactions').

''Further discussion'': 
What effects have been neglected in these calculations of Diels Alder transition states?

Look at published examples and investigate further if you have time.
(e.g. {{DOI|10.1021/jo0348827}})

----
See also: [[Mod:timetable|Timetable]],[[Mod:lectures|Intro lecture]], [[mod:programs|Programs]], [[mod:organic|Module 1]], [[Mod:inorganic|Module 2]], [[Mod:phys3|Module 3]]

© 2008-2011, Imperial College London

Gaussian Bugs and Issues

2013-01-09T14:45:47Z

Lmt09: /* Modredundant */

==Gaussian Bugs and Issues==

Section to place any suspected bugs and errors

===CASSCF===

'''1.'''Quartic linear search updated in l103 seems to change which inital CASSCF guess is read (Lee).

Use IOp(4/5=1) to read from the checkpoint file. 
Only a problem with GDV versions (GDVH13 and GDVH23 checked), not G09 (G09A02 and G09C01 checked)

'''2.'''Freq and StateAverage kewords are incompatible. Fix needs to use Freq and IOp for SA (David).

===CIS===

'''1.'''Geomtry optimization gives a core dump when reading initial guess vectors from the rwf. (Lee)

Use IOp(9/49=3) to force read from the checkpoint file. 
The checkpoint file must be available to read from in the first step 
If this is not the case run a single point and then run the optimization using the SP chk file with the IOp 
Checked with GDVH13 only.

===Opt===

====Modredundant====

'''1.''' Using the Frozen coordinate command whilst setting the value for the coordinate (i.e. B 1 2 2.2 F which should freeze the bond distance between atoms 1 and 2 to 2.2 angstrom) sets the correct angle but does not freeze it. (Lee)

Build the coordinate and freeze in different specifications e.g. 
B 1 2 2.2 B 
B 1 2 F

Gaussian Bugs and Issues

2013-01-09T14:45:29Z

Lmt09:

==Gaussian Bugs and Issues==

Section to place any suspected bugs and errors

===CASSCF===

'''1.'''Quartic linear search updated in l103 seems to change which inital CASSCF guess is read (Lee).

Use IOp(4/5=1) to read from the checkpoint file. 
Only a problem with GDV versions (GDVH13 and GDVH23 checked), not G09 (G09A02 and G09C01 checked)

'''2.'''Freq and StateAverage kewords are incompatible. Fix needs to use Freq and IOp for SA (David).

===CIS===

'''1.'''Geomtry optimization gives a core dump when reading initial guess vectors from the rwf. (Lee)

Use IOp(9/49=3) to force read from the checkpoint file. 
The checkpoint file must be available to read from in the first step 
If this is not the case run a single point and then run the optimization using the SP chk file with the IOp 
Checked with GDVH13 only.

===Opt===

====Modredundant====

'''1.''' Using the Frozen coordinate command whilst setting the value for the coordinate (i.e. B 1 2 2.2 F which should freeze the bond distance between atoms 1 and 2 to 2.2 angstrom) sets the correct angle but does not freeze it. (Lee)

Build the coordinate and freeze in different specifications e.g. 
B 1 2 2.2 B
B 1 2 F

Gaussian Bugs and Issues

2013-01-09T12:34:13Z

Lmt09: /* CASSCF */

Gaussian Bugs and Issues

2012-12-10T15:32:29Z

Lmt09: /* CASSCF */

Gaussian Bugs and Issues

2012-11-29T12:20:01Z

Lmt09:

==Gaussian Bugs and Issues==

Section to place any suspected bugs and errors

===CASSCF===

Quartic linear search updated in l103 seems to change which inital CASSCF guess is read (Lee).

Freq and StateAverage kewords are incompatible. Fix needs to use Freq and IOp for SA (David).

Gaussian Bugs and Issues

2012-11-29T12:19:41Z

Lmt09: Created page with "==Gaussian Bugs and Issues== Section to place any suspected bugs and errors ===CASSCF=== Quartic linear search updated in l103 seems to change which inital CASSCF guess is rea..."

==Gaussian Bugs and Issues==

Section to place any suspected bugs and errors

===CASSCF===

Quartic linear search updated in l103 seems to change which inital CASSCF guess is read (Lee).
Freq and StateAverage kewords are incompatible. Fix needs to use Freq and IOp for SA (David).

Resgrp:comp-photo

2012-11-29T12:10:47Z

Lmt09: /* Example Calculations / Tutorials */

= Computational Photochemistry Research Group Wiki, Imperial College London =

A place to write something about computing facilities (e.g. CX1),
documenting codes (Gaussian + dynamics; latest versions)
and example input files.
-[[User:Mjbear|Mjbear]] 16:16, 17 October 2008 (BST)

Please amend and update this! Log in using your IC account.

== Computing Resources ==
=== [[Resgrp:comp-photo-hpc|Using College HPC to run Gaussian: CX1]] ===
===== [[Resgrp:comp-photo-hpc-ax1|- AX1 / AX2 ]] =====
=== [[Resgrp:utilities|Utilities]] ===

== Codes ==
=== [[Resgrp:comp-photo-gaussian|Gaussian versions]] ===
===== [[Resgrp:comp-photo-gaussian_problems|- known problems ]] =====
===Developing Gaussian===
There's a new regime!
====[[Resgrp:comp-photo-version control-local-dev-master|Local development versions]]====
A page explaining what is in the local development repository. Currently rather sparse,
you are encouraged to be much more verbose.
====[[Resgrp:comp-photo-version control|Using Mercurial]]====
A rather long manual on why and how to use the mercurial version control system,
somewhat tailored to Gaussian development.
====[[Resgrp:comp-photo-version control short|Quick guide to Mercurial]]====
A short version of the above.
====[[Resgrp:comp-photo-version control example|Example workflows]]====
A couple of examples of how one would use mercurial to develop Gaussian
====[[Resgrp:comp-photo-new gdv layout|The new layout]]====
An explanation of the new layout of the Gaussian source code necessitated
by version control. Plus some tips and tricks.
=== [[Resgrp:comp-photo-dyn|MCTDH (DD-vMCG)]] ===
=== [[Resgrp:comp-photo-onthefly|Direct CAS/RAS]] ===

== Example Calculations / Tutorials ==
=== [[Resgrp:comp-photo-benzene-tutorial|Benzene CASSCF tutorial (G03)]] ===
=== [[Resgrp:comp-photo-CHD|Seam calculations of CHD]] ===
===== [[Resgrp:comp-photo-seam|- Current Seam Frequency Code ]] =====
=== [[Resgrp:comp-photo/sh|Classical path trajectories]] ===

=== [[Resgrp:comp-photo/VB|Valence-Bond analysis]] ===

=== [[Resgrp:comp-photo-visualise|GaussView, visualisation etc]] ===

=== [[ONIOM]] ===

=== [[Gaussian Bugs and Issues]] ===

== Calculations ==
=== [[Resgrp:comp-photo-calculations|Examples]] ===

Resgrg:comp-photo-version control-codes

2012-11-26T16:46:15Z

Lmt09: /* ONIOM Freq */

This page explains the local modifications made to the Gaussian code by the group. Hopefully it'll also explain
how to use them as well.

==The codes==

===Shaopeng===
This is a replacement for the existing matrix element generator (hereafter known as the Klene code)
that is used for large (number of orbitals > 8)
CASSCF calculations. If your job shows <code>SME calculated on fly</code> then it's using this path.

The Shaopeng code computes the same matrix elements but does it faster, and in a way that allows certain compiler
optimisations to be used. The result should be an approximate 4X speedup, depending on the job and the architecture that you're running on. This speedup comes at a cost of memory. The algorithm is much more demanding of memory and this extra demand scales with the number of processors.

How to use it: by default l510 will make a decision on whether to use the Klene code path or the Shaopeng based on the problem size and the amount of memory and processors allocated. You can force one code path over the other
with the IOP(5/139=x) where:

{|
! Iop(5/139)
! Effect
|-
| 0
| Default, choose code path automatically
|-
| 1
| Force Klene code path
|-
| 2
| Force Shaopeng code path
|}

===Seamfollowing===

Dunno

===ONIOM Freq===

Shows contributions of sub-calculations to the ONIOM frequency according to the formula:

<math>(f^{ONIOM}_{i})^{2} = (f^{'high}_{model,i})^{2} + (f^{'low}_{real,i})^{2}-(f^{'low}_{model,i})^{2}</math>

The sub-calculation terms are obtained from the diagonal elements after transformation of the sub-Hessians with the ONIOM normal modes.
This gives the curvature of each sub-calculation PES with respect to the ONIOM mode so we can see how the ONIOM PES is determined
i.e. if a transition state located using ONIOM has an imaginary mode comprised of three real sub-components, then the transition state is a result of the extrapolation and is not consistent with the underlying calculations. This result suggests caution should be used as the stationary point may be spurious.

The ONIOM frequency analysis is called using IOp(7/123=1):

IOp(123) ... Print partitioning of ONIOM vibrational frequencies
into contributions from individual sub-calculations.
see Vreven et al. JCTC, 2012, DOI: 10.1021/ct300612m
0 ... Don't do ONIOM frequency analysis
1 ... Do ONIOM frequency analysis

Example output:

Mode Low-Real Low-Model High-Model Total
1 2181.3672 433.6860 231.1949 2150.2862
2 4029.7184 2543.4278 2310.9318 3887.1598
3 4317.3046 869.1994 789.5792 4301.9818

Mod:phys3

2012-10-29T18:03:37Z

Lmt09: /* Optimizing the Reactants and Products */

See also:[[Mod:timetable|Timetable]], [[Mod:lectures|Intro lecture]], [[mod:programs|Programs]], [[mod:organic|Module 1]], [[Mod:inorganic|Module 2]], [[Mod:phys3|Module 3]], [http://www.gaussian.com/g_tech/gv5ref/gv5ref_toc.htm Gaussian Online User Manual] | [http://faculty.ycp.edu/~jforesma/educ/visual/index.html Visualization Tutorials]
= Module 3 =

In this set of computational experiments, you will characterise transition structures on potential energy surfaces for the Cope rearrangement and Diels Alder cycloaddition reactions.

There are two parts:
a) tutorial material: how to use the programs and methods,
b) more challenging examples, with guidelines but fewer explicit instructions.



In the second year physical chemistry laboratory, you may have carried out dynamics calculations using model potential energy surfaces to explore transition states. In that computational experiment, the total energy could quickly be calculated for different geometries of a triatomic system using an analytical function of the atomic coordinates (for more information, see for example [http://books.google.com/books?id=T8IZ1aa_FRkC&pg=RA1-PA36&lpg=RA1-PA36&dq=%22lake+eyring%22&source=web&ots=OXY00lSZ7D&sig=Ld_MTNwNjUDNGzB_5w1IxaMBMPU&hl=en&sa=X&oi=book_result&resnum=7&ct=result here] and [http://www.rsc.org/ejarchive/DC/1979/DC9796700007.pdf here]).

In this experiment, you will be studying transition structures in larger molecules. There are no longer fitted formulae for the energy, and the molecular mechanics / force field methods that work well for structure determination cannot be used (in general) as they do not describe bonds being made and broken, and changes in bonding type / electron distribution. (This is the main difference from Module 1). Instead, we use molecular orbital-based methods, numerically solving the Schrodinger equation, and locating transition structures based on the local shape of a potential energy surface. As well as showing what transition structures look like, reaction paths and barrier heights can also be calculated.

==The Cope Rearrangement Tutorial==



''This part of the module is described as a 'tutorial' because it's an introduction to various computational techniques for locating transition structures on potential energy surfaces. It's different to the GaussView tutorial you may have worked through earlier: it's an exercise where you're given specific instructions, see if you can follow them, and also whether there are problems or better ways of carrying the exercise out. Please include this part in your write-up. Marks will be given for correct answers, the documentation showing how you got these, discussion, and how you went about solving any problems you encountered.''

In this tutorial we will use the Cope rearrangement of 1,5-hexadiene as an example of how to study a chemical reactivity problem.

Your objectives are to locate the low-energy minima and transition structures on the C6H10 potential energy surface, to determine the preferred reaction mechanism.

[[Image:pic1.jpg|right|thumb|Cope rearrangement]]

This [3,3]-sigmatropic shift rearrangement has been the subject of numerous experimental and computational studies (e.g. Houk et al. {{DOI|10.1021/ja00101a078}}), and for a long time its mechanism (concerted, stepwise or dissociative) was the subject of some controversy. Nowadays it is generally accepted that the reaction occurs in a concerted fashion via either a "chair" or a "boat" transition structure, with the "boat" transition structure lying several kcal/mol higher in energy. The B3LYP/6-31G* level of theory has been shown to give activation energies and enthalpies in remarkably good agreement with experiment. In this tutorial we will show how these can be calculated using Gaussian.

 

{| align="center"
|- align="center"
|
[[Image:pic2a.jpg]]
|
[[Image:pic2b.jpg]]
|- align="center"
| align="center" | ''Chair Transition State''
| align="center" | ''Boat Transition State''
|}

 

===Optimizing the Reactants and Products===

'' In this section you will learn how to optimize a structure, symmetrize it to find its point group, calculate and visualize vibrational frequencies and correct potential energies in order to compare them with experimental values. It is assumed that you are already familiar with using the builder in GaussView. ''

 

(a) Using GaussView, draw a molecule of 1,5-hexadiene with an "anti" linkage (aproximately a.p.p conformation) for the central four C atoms . Clean the structure using the '''Clean''' function under the '''Edit''' menu.

Now we will optimize the structure at the '''HF/3-21G''' level of theory. Select '''Gaussian''' under the '''Calculate''' menu, click on the '''Job Type''' tab and choose '''Optimization'''. The default method should already be Hartree Fock and the default basis set is 3-21G, so there should be no need to change these. You can check this by clicking on the '''Method''' tab. Change the %mem under the '''Link 0''' tab to 250 MB (though this could be increased to 500 MB). Submit the job by clicking on the '''Submit''' button at the bottom of the window and give the job a meaningful name (e.g. react_anti).

When the job has finished, you will be asked if you want to open a file. Select '''Yes''' and choose the checkpoint (chk) file with the name of the job you have just run (e.g. react_anti.chk). This checkpoint file is a binary file that stores data calculated by Gaussian. The name of the chk file should have been assigned by default, but by default, this file will be created in the C:\Windows\G03\Scratch folder. Once the file has been opened, click on the '''Summary''' button under the '''Results''' menu and make a note of the energy.


Does your final structure have symmetry? Select '''Symmetrize''' under the '''Edit''' menu (note that sometimes it is necessary to relax the search criteria under the '''Point Group''' menu). Make a note of the point group.

 

(b) Now draw another molecule of 1,5-hexadiene with a "gauche" linkage for the central four C atoms. Would you expect this structure to have a lower or a higher energy than the anti structure you have just optimized? Optimize the structure at the '''HF/3-21G''' level of theory and compare your final energy with that obtained in (a). Again, check if the molecule has symmetry and make a note of the point group.

 

(c) Normally, calculated activation energies and enthalpies use the lowest energy conformation of a reactant molecule as a reference. Based on your results from above, try to predict what the lowest energy conformation of 1,5-hexadiene might be. Test out your hypothesis by drawing the structure and optimizing it.

 

(d) A table containing the low energy conformers of 1,5-hexadiene and their point groups is shown in [[Mod:phys3#Appendix 1|Appendix 1]]. Compare the structures that you have optimized with those in the table and see if you can identify your structure.

 

(e) Draw the Ci ''anti2'' conformation of 1,5-hexadiene (unless you have already located it). Optimize it at the '''HF/3-21G''' level of theory and make sure it has Ci symmetry. Compare your final energy to the one given in the table.


 

(f) When you are happy that your structure is the same as the one in the table, reoptimize it at the '''B3LYP/6-31G*''' level (6-31G* is equivalent to 6-31G(d) by selecting '''DFT''' under the '''Method''' menu and '''B3LYP''' from the box with the functionals on the right-hand side. Now select '''Link 0''' and change the name of the chk file to the name of the DFT optimization that you are about to run. Note that it is always advisable to do this when re-using or modifying existing structures to ensure that the original chk file is not overwritten. Run the job and make a note of the energy. Now compare the final structures from the '''HF/3-21G''' calculation with that at the higher level of theory. How much does the overall geometry change?

 

(g) The final energies given in the output file represent the energy of the molecule on the bare potential energy surface. To be able to compare these energies with experimentally measured quantities, they need to include some additional terms, which requires a frequency calculation to be carried out. The frequency calculation can also be used to characterize the critical point, i.e. to confirm that it is a minimum in this case: that all vibrational frequencies are real and positive.

Starting from your optimized B3LYP/6-31G* structure, run a frequency calculation at the same level of theory. You can do this by selecting '''Frequency''' under the '''Job Type''' tab. Ensure that the method is still correctly specified under the '''Method''' tab (''caution: on Windows, sometimes 'scrf=(solvent=water,check)' is incorrectly added!'') and then change the name of the chk file under the '''Link 0''' tab to the name of the frequency job that you are about to run. Run the job. Once the job has finished, open the log file this time. Select '''Vibrations''' under the '''Results''' menu. A list of all the vibrational frequencies modes should appear. Check that there are no imaginary frequencies, only real ones. You can visualize some of these vibrations under this menu and simulate the infrared spectrum.


Now, select '''View File''' under the '''Results''' menu and open the output file in the visualizer. Scroll down to the section beginning '''Thermochemistry'''. Under the vibrational temperatures a list of energies should be printed. Make a note of (i) the sum of electronic and zero-point energies, (ii) the sum of electronic and thermal energies, (iii) the sum of electronic and thermal enthalpies, and (iv) the sum of electronic and thermal free energies. The first of these is the potential energy at 0 K including the zero-point vibrational energy (E = Eelec + ZPE), the second is the energy at 298.15 K and 1 atm of pressure which includes contributions from the translational, rotational, and vibrational energy modes at this temperature (E = E + Evib + Erot + Etrans), the third contains an additional correction for RT (H = E + RT) which is particularly important when looking at dissociation reactions, and the last includes the entropic contribution to the free energy (G = H - TS). It is important to make sure that you select the correct energy/enthalpy term to compare to your experimental values. Note that these corrections can also be calculated at other temperatures using the '''Temperature''' option in Gaussian, If you have time, try re-calculate these quantities at 0 K as shown in the [[mod:gv_advanced | Advanced GaussView Tutorial]].

 

===Optimizing the "Chair" and "Boat" Transition Structures ===

'' In this section you will learn how to set up a transition structure optimization (i) by computing the force constants at the beginning of the calculation, (ii) using the redundant coordinate editor, and (iii) using QST2. You will also visualize the reaction coordinate and run the IRC (Intrinisic Reaction Coordinate) and calculate the activation energies for the Cope rearrangement via the "chair" and "boat" transition structures. ''

 

The "chair" and "boat" transition structures for the Cope rearrangement are shown in [[Mod:phys3#Appendix 2|Appendix 2]]. Both consist of two C3H5 allyl fragments positioned approximately 2.2 Å apart, one with C2h symmetry and the other with C2v symmetry.

 

(a) Draw an allyl fragment (CH2CHCH2) and optimize it using the '''HF/3-21G''' level of theory. Your structure should look like one half of the transition structures shown below.

Now open a new GaussView window by going to the '''File''' menu and selecting '''New''' and then '''Create MolGroup'''. Copy the optimized allyl structure from the first calculation by selecting '''Copy''' under the '''Edit''' menu, and then paste it twice into the new window by selecting '''Paste''' and then '''Append Molecule'''. Now orient the two fragments so that they look roughly like the chair transition state below by using the '''Shift Alt keys + Left Mouse button''' to translate one fragment with respect to the other and the '''Alt key + Left Mouse button''' to rotate it. The distance between the terminal ends of the allyl fragments should be approximately 2.2 Å apart. Save this structure to a Gaussian input file with a meaningful name (e.g. chair_ts_guess).

We are now going to optimize this transition state manually in two different ways. Transition state optimizations are more difficult than minimizations because the calculation needs to know where the negative direction of curvature (i.e. the reaction coordinate) is. If you have a reasonable guess for your transition structure geometry, then normally the easiest way to produce this information is to compute the force constant matrix (also known as the Hessian) in the first step of the optimization which will then be updated as the optimization proceeds. This is what we will try to do in the next section. However, if the guess structure for the transition structure is far from the exact structure, then this approach may not work as the curvature of the surface may be significantly different at points far removed from the transition structure. In some cases, a better transition structure can be generated by freezing the reaction coordinate (using '''Opt=ModRedundant''' and minimizing the rest of the molecule. Once the molecule is fully relaxed, the reaction coordinate can then be unfrozen and the transition state optimization is started again. One advantage of doing this, is that it may not be necessary to compute the whole Hessian once this has been done, and just differentiating along the reaction coordinate might give a good enough guess for the initial force constant matrix. This can save a considerable amount of time in cases where the force constant calculation is expensive.

 

(b) Use Hartree Fock and the default basis set 3-21G for parts (b) to (f).

Create a new MolGroup ('''File → New → Create MolGroup''') and copy and paste your guess structure into the window. Now set up a Gaussian optimization for a transition state. Go to the '''Gaussian''' menu under '''Calculate''' and click on the '''Job Type''' tab. Select '''Opt+Freq''' and then change '''Optimization to a Minimum''' to '''Optimization to a TS (Berny)'''. Choose to calculate the force constants '''Once''' and in the Additional keyword box at the bottom, type '''Opt=NoEigen'''. The latter stops the calculation crashing if more than one imaginary frequency is detected during the optimization which can often happen if the guess transition structure is not good enough. Submit the job. If the job completes successfully, you should have optimized to the structure shown in [[Mod:phys3#Appendix 2|Appendix 2]] and the frequency calculation should give an imaginary frequency of magnitude 818 cm-1. Animate the vibration and ensure that it is the one corresponding to the Cope rearrangement.

 

(c) Now we will try optimizing the transition structure again using the frozen coordinate method. Create a new MolGroup ('''File → New → Create MolGroup''') and copy and paste your guess structure into the window again. Now select '''Redundant Coord Editor''' from the '''Edit''' menu. Click on the highlighted file icon at the top left-hand corner (Create a New Coordinate) and a line should appear below saying '''Add Unidentified (?, ?, ?, ?)'''. Now go back to the GaussView window and select two of the terminal carbons from the allyl fragments which form/break a bond during the rearrangement. Return to the coordinate editor and select '''Bond''' instead of '''Unidentified''' and select '''Freeze Coordinate''' instead of '''Add'''. Now click on the icon again to generate another coordinate. This time select the opposite two terminal atoms and again select '''Bond''' and '''Freeze Coordinate'''. Click OK. Now set up the optimization as if it were a minimum and you should see the option '''Opt=ModRedundant''' already included in the input line. Submit the job.

'''Note:''' GaussView allows you to produce an input file with the frozen coordinate specified as e.g. <tt>B 5 1 2.200000 F</tt>. Unfortunately, a recent update to the Gaussian program means it does not recognise this syntax, and just ignores this line. This means that the coordinate ends up being optimised rather than frozen. Therefore do not use this method, but ensure the guess structure has suitable guess transition bond distances(~2.2 Å) using the ''Modify Bond'' tool in GaussView --[[User:Rzepa|Rzepa]] 14:39, 29 October 2012 (UTC)

 

(d) When the job has finished, open the chk file. You should find that the optimized structure looks a lot like the transition you optimized in section (b), except the bond forming/breaking distances are fixed to 2.2 Å. Now we are going to optimize them too. Open the '''Redundant Coord Editor''' from the '''Edit''' menu again and create a new coordinate as before by clicking on the icon, Select one of the bonds that was previously frozen and this time choose '''Bond''' instead of '''Unidentified''' and '''Derivative''' instead of '''Add'''. Repeat the procedure for the other bond. This time you need to set up a transition state optimization but we are not going to calculate the force constants as we did in section (b) (so we leave this option as '''Never'''), instead we will use a normal guess Hessian modified to include the information about the two coordinates we are differentiating along. Change the name of the chk file in '''Link 0''' if you do not want to write over the previous calculation and submit the job. When the calculation has finished, open the chk file, check the bond forming/bond breaking bond lengths and compare the structure to the one you optimized in section (b).

 

(e) Now we will optimize the boat transition structure. We will do this using the '''QST2''' method. In this method, you can specify the reactants and products for a reaction and the calculation will interpolate between the two structures to try to find the transition state between them. You must make sure that your reactants and products are numbered in the same way. Therefore, although our reactants and products are both 1,5-hexadiene, we will need to manually change the numbering for the product molecule so that it corresponds to the numbering obtained if our reactant had rearranged.

e.g.

<center>[[Image:pic3.jpg|200px]]</center>

Open the chk file corresponding to the optimized Ci reactant molecule (''anti2'' in [[Mod:phys3#Appendix 1|Appendix 1]]). Now open a second window and create a new MolGroup. Copy the optimized reactant molecule into the new window. In the same window, now select '''File → New → Add to MolGroup'''. The original molecule should disappear and a green circle should appear at the top left-hand corner with a '''2''' next to it. Clicking on the down arrow by the '''2''' will take you back to the original window and you will see your molecule again. This is how we read multiple geometries into GaussView. Go back to window '''2''', and copy and paste the reactant molecule a second time. This is going to be the product molecule and will be the molecule on which we need to change the numbering. If you now click on the icon showing two molecules side by side, then you can view both molecules simultaneously.

Now go to the '''View''' menu and select '''Labels''' so that you can see the numbering on both structures. Orient the two structures separately so they look something like the following:

{| align="center"
|- align="center"
|
[[Image:pic4a.jpg|200px]]
|
[[Image:pic4b.jpg|200px]]
|- align="center"
| align="center" | ''Reactant''
| align="center" | ''Product''
|}

 

Now click on the product structure. Go to the '''Edit''' menu and select '''Atom List'''. Starting from Atom 1 on the reactant, go through and renumber all the atoms on the Product so that they match the reactant molecule, e.g. for the numbering above you would start by changing atom '''6''' on the product molecule to atom '''1'''. The other atom numbers will update as you do this so make sure you do it in the correct order. At the end, the numbering on your two molecules should correspond to each other in the following way:

 

{| align="center"
|- align="center"
|
[[Image:pic5a.jpg|200px]]
|
[[Image:pic5b.jpg|200px]]
|- align="center"
| align="center" | ''Reactant''
| align="center" | ''Product''
|}

 

Now we will set up the first '''QST2''' calculation. Go to the '''Gaussian''' menu and select '''Job Type''' as '''Opt+Freq''', and optimize to a transition state. This time you will have two options - '''TS (Berny)''' which we used in the previous calculations and '''TS (QST2)'''. Select '''TS (QST2)'''. Submit the job.

You will find that the job fails. To see why, open the chk file you created and view the structure. You will see that it looks a bit like the chair transition structure but more dissociated. In fact when the calculation linearly interpolated between the two structures, it simply translated the top '''allyl''' fragment and did not even consider the possibility of a rotation around the central bonds. It is clear that the QST2 method is never going to locate the boat transition structure if we start from these reactant and product structures.

Now go back to the original input file where you set up your QST2 calculation. We will now modify the reactant and product geometries so that they are closer to the boat transition structure. Click on the reactant molecule first and select the central '''C-C-C-C''' dihedral angle (i.e. '''C2-C3-C4-C5''' for the molecule above) and change the angle to 0o. Then select the inside '''C-C-C''' (i.e. '''C2-C3-C4''' and '''C3-C4-C5''' for the molecule above) and reduce them to 100o. Do the same for the product molecule. Your reactant and product molecules should now look like the following:

 

{| align="center"
|- align="center"
|
[[Image:pic6a.jpg|200px]]
|
[[Image:pic6b.jpg|200px]]
|- align="center"
| align="center" | ''Reactant''
| align="center" | ''Product''
|}

 

Set up the QST2 calculation again, renaming both the chk file under '''Link 0''' and the input file. Run the job again. This time it should converge to the boat transition structure. Check that there is only one imaginary frequency and visualize its motion.

The object of this exercise is to illustrate that although the QST2 method is has some advantages because it is fully automated, it can often fail if your reactants and products are not close to the transition structure. There is another method, the '''QST3''' method, that allows you to input the geometry of a guess transition structure also and this can often be more reliable. If you have time, you can try generating a guess boat transition structure and see if you can get the calculation to converge using the original reactant and product molecules. Remember to check the atom numbers in the transition structure are in the right order.

 

(f) Take a look at your optimized chair and boat transition structures. Which conformers of 1,5-hexadiene do you think they connect? You will find that it is almost impossible to predict which conformer the reaction paths from the transitions structures will lead to. However, there is a method implemented in Gaussian which allows you to follow the minimum energy path from a transition structure down to its local minimum on a potential energy surface. This is called the '''Intrinsic Reaction Coordinate''' or '''IRC''' method. This creates a series of points by taking small geometry steps in the direction where the gradient or slope of the energy surface is steepest.

Open the chk file for one of your optimized chair transition structures. Under the '''Gaussian''' menu, select '''IRC''' under the '''Job Type''' tab. You will be presented with a number of options. The first is to decide whether to compute the reaction coordinate in one or both directions. As our reaction coordinate is symmetrical, we will only choose to compute it in the forward direction. Normally you would do both forward and reverse, either in one job or in two separate jobs. You are also given the option to calculate the force constants once, at every step along the IRC or to read them from the chk file. You would use the latter option if you have previously run a frequency calculation. In this case, to avoid confusion with chk files, we will just recompute them at the beginning of the calculation. (The '''IRCMax''' option can also be specified here. This takes a transition structure as its input, and finds the maximum energy along a specified reaction path, taking into account zero-point energy etc., and produces all the quantities needed for a variational transition state theory calculation. We will leave this unchecked for the purposes of this exercise.) The final option to consider is the number of points along the IRC. The default is '''6''' but this is normally never enough. Let's change this to 50 and see how the calculation progresses. Change the name of the chk file under '''Link 0''' and submit the job. The job will take a while so now is a good time to take a coffee break...

When the IRC calculation has finished, open the chk file with all the intermediate geometries and see how the calculation has progressed. You will find that it hasn't reached a minimum geometry yet. This leaves you three options: (i) you can take the last point on the IRC and run a normal minimization; (ii) you can restart the IRC and specify a larger number of points until it reaches a minimum; (iii) you can redo the IRC specifying that you want to compute the force constants at every step. There are advantages and disadvantages to each of these approaches. Approach (i) is the fastest, but if you are not close enough to a local minimum, you may end up in the wrong minimum. Approach (ii) is more reliable but if too many points are needed, then you can also veer off in the wrong direction after a while and end up at the wrong structure. Approach (iii) is the most reliable but also the most expensive and is not always feasible for large systems. You can try any or all of these approaches and see which conformation you end up in.

 

(g) Finally we need to calculate the activation energies for our reaction via both transition structures. To do this we will need to reoptimize the chair and boat transition structures using the '''B3LYP/6-31G*''' level of theory and to carry out frequency calculations. You can start from the HF/3-21G optimized structures. Once the calculations have converged, compare both the geometries and the difference in energies between the reactants and transition states at the two levels of theory. What you should find is that the geometries are reasonably similar, but the energy differences are markedly different.

As a consequence of this, it is often more computational efficient to map the potential energy surface using the low level of theory first and then to reoptimize at the higher level as we have done in this exercise.

The experimental activation energies are 33.5 ± 0.5 kcal/mol via the chair transition structure and 44.7 ± 2.0 kcal/mol via the boat transition structure at 0 K. If you take the values computed at 0 K, how close are they to the experimental values? You can also find the energies with thermal correction at 298.15 K under the Thermochemistry data in the output file. If you have time, you can recompute them at higher temperature. Alternatively, you can use the utility program '''FreqChk''' to obtain energies at a different temperature. This only requires the chk file from a frequency calculation and allows you to retrieve frequency and thermochemistry data as well as calculating them with an alternate temperature, pressure, scale factor, and/or isotope substitutions. The '''FreqChk''' utility program can be accessed from '''Gaussian03W'''. Launch '''Gaussian03W'''. Select '''utilities''' from the menu and click on '''FreqChk''' to launch the utility program. You will be prompted for a chk file. Select your chk file from the C:\G03W\Scratch directory and follow the instructions from this [http://www.gaussian.com/g_tech/g_ur/u_freqchk.htm web link] to proceed.

 

=== Appendix 1 ===

 

{| border="1" cellpadding="5"
|- align="center"
| width="150" | '''Conformer'''
| width="150" | '''Structure'''
| width="100" | '''Point Group'''
| width="200" | '''Energy/Hartrees HF/3-21G'''
| width="200" | '''Relative Energy/kcal/mol'''
|- align="center"
| ''gauche1''
|
[[Image:gauche1.jpg|150px]]
| C2
| -231.68772
| 3.10
|- align="center"
| ''gauche2''
|
[[Image:gauche2.jpg|150px]]
| C2
| -231.69167
| 0.62
|- align="center"
| ''gauche3''
|
[[Image:gauche3.jpg|150px]]
| C1
| -231.69266
| 0.00
|- align="center"
| ''gauche4''
|
[[Image:gauche4.jpg|150px]]
| C2
| -231.69153
| 0.71
|- align="center"
| ''gauche5''
|
[[Image:gauche5.jpg|150px]]
| C1
| -231.68962
| 1.91
|- align="center"
| ''gauche6''
|
[[Image:gauche6.jpg|150px]]
| C1
| -231.68916
| 2.20
|- align="center"
| ''anti1''
|
[[Image:anti1.jpg|150px]]
| C2
| -231.69260
| 0.04
|- align="center"
| ''anti2''
|
[[Image:anti2.jpg|150px]]
| Ci
| -231.69254
| 0.08
|- align="center"
| ''anti3''
|
[[Image:anti3.jpg|150px]]
| C2h
| -231.68907
| 2.25
|- align="center"
| ''anti4''
|
[[Image:anti4.jpg|150px]]
| C1
| -231.69097
| 1.06
|}

 

=== Appendix 2 ===

{| cellpadding="5"
|- align="center"
|
[[Image:appendix2a.jpg|300px]]
|- align="center"
| ''C2h Chair Transition State''
|- align="center"
|- align="center"
|- align="center"
|
[[Image:appendix2b.jpg|300px]]
|- align="center"
| ''C2v Boat Transition State''
|}

 

=== Results Table ===

 

''' Summary of energies (in hartree) '''

{| border="1" cellspacing="1" cellpadding="10"
|-
| ''' '''
!colspan="3"|'''HF/3-21G'''
!colspan="3"|'''B3LYP/6-31G*'''
|-
| ''' '''
| width="125" align="center" | '''Electronic energy'''
| width="125" align="center" | '''Sum of electronic and zero-point energies'''
| width="125" align="center" | '''Sum of electronic and thermal energies'''
| width="125" align="center" | '''Electronic energy'''
| width="125" align="center" | '''Sum of electronic and zero-point energies'''
| width="125" align="center" | '''Sum of electronic and thermal energies'''
|-
| ''' '''
| ''' '''
| width="125" align="center" | '''at 0 K'''
| width="125" align="center" | '''at 298.15 K'''
| ''' '''
| width="125" align="center" | '''at 0 K'''
| width="125" align="center" | '''at 298.15 K'''
|-
| '''Chair TS'''
| align="center" | -231.619322
| align="center" | -231.466705
| align="center" | -231.461346
| align="center" | -234.556983
| align="center" | -234.414919
| align="center" | -234.408998
|-
| '''Boat TS'''
| align="center" | -231.602802
| align="center" | -231.450929
| align="center" | -231.445300
| align="center" | -234.543093
| align="center" | -234.402340
| align="center" | -234.396006
|-
| '''Reactant (''anti2'')'''
| align="center" | -231.692535
| align="center" | -231.539539
| align="center" | -231.532566
| align="center" | -234.611710
| align="center" | -234.469203
| align="center" | -234.461856
|}

 *1 hartree = 627.509 kcal/mol 

''' Summary of activation energies (in kcal/mol) '''

{| border="1" cellspacing="1" cellpadding="10"
|-
| ''' '''
| width="125" align="center" | '''HF/3-21G'''
| width="125" align="center" | '''HF/3-21G'''
| width="125" align="center" | '''B3LYP/6-31G*'''
| width="125" align="center" | '''B3LYP/6-31G*'''
| width="125" align="center" | '''Expt.'''
|-
| ''' '''
| width="125" align="center" | '''at 0 K '''
| width="125" align="center" | '''at 298.15 K'''
| width="125" align="center" | '''at 0 K'''
| width="125" align="center" | '''at 298.15 K'''
| width="125" align="center" | '''at 0 K'''
|-
| '''ΔE (Chair)'''
| align="center" | 45.70
| align="center" | 44.69
| align="center" | 34.06
| align="center" | 33.17
| align="center" | 33.5 ± 0.5
|-
| '''ΔE (Boat)'''
| align="center" | 55.60
| align="center" | 54.76
| align="center" | 41.96
| align="center" | 41.32
| align="center" | 44.7 ± 2.0
|}

 

==The Diels Alder Cycloaddition==

In this exercise, you will characterise transition structures using any of the methods described above in the tutorial: the choice is up to you. In addition, you will look at the shape of some of the molecular orbitals. To help you structure your report, there is a data/discussion sheet at the end of this section.

[[Image:mb_da1.jpg |right|thumb|Diels Alder cycloaddition]]
The Diels Alder reaction belongs to a class of reactions known as pericyclic reactions. The π orbitals of the dieneophile are used to form new σ bonds with the π orbitals of the diene. Whether or not the reactions occur in a concerted stereospecific fashion ('''allowed''') or not ('''forbidden''') depends on the number of π electrons involved. In general the HOMO/LUMO of one fragment interacts with the HOMO/LUMO of the other reactant to form two new bonding and anti-bonding MOs. The nodal properties allow one to make predictions according to the following rule:

''If the HOMO of one reactant can interact with the LUMO of the other reactant then the reaction is '''allowed'''.''

''The HOMO-LUMO can only interact when there is a significant overlap density. If the orbitals have different symmetry properties then no overlap density is possible and the reaction is '''forbidden'''.''

If the dieneophile is substituted, with substituents that have π orbitals that can interact with the new double bond that is being formed in the product, then this interaction can stabilise the regiochemistry (i.e. head to tail versus tail to head) of the reaction. In this exercise you will study the nature of the transition structure for the Diels Alder reaction, both for the prototypical reaction and for the case where both diene and dieneophile carry substituents, and where secondary orbital effects are possible. Clearly, the factors that control the nature of the transition state are quantum mechanical in origin and thus we shall use methods based upon quantum chemistry.

 

Shown on the right is a diagram of the transition state for the Diels-Alder reaction between ethylene and butadiene. The ethylene approaches the cis form of butadiene from above.
[[Image:mb_da2.jpg |right|thumb|Ethylene+Butadiene cycloaddition]]

Before beginning our quantitative study, it is helpful to discuss the interaction of the π orbitals in a simple qualitative way.

''You will confirm some of these considerations in your computations.''

The principal orbital interactions involve the π/ π* orbitals of ethylene and the HOMO/LUMO of butadiene. It is referred to as [4s + 2s] since one has 4 π orbitals in the π system of butadiene. The orbitals of ethylene and butadiene and ethylene can be classified as symmetric '''s''' or anti-symmetric '''a''' with respect to the plane of symmetry shown.

The HOMO of ethylene and the LUMO of butadiene are both '''s''' (symmetric with respect to the reflection plane) and the LUMO of ethylene and the HOMO of butadiene are both '''a'''. Thus it is the HOMO-LUMO pairs of orbital that interact, and energetically, the HOMO of the resulting adduct with two new σ bonds is '''a'''.

 

=== Exercise ===

Use the the AM1 semi-empirical molecular orbital method for these calculations (to start with).

i) Use GaussView to build cis butadiene, and optimize the geometry using Gaussian. Plot the HOMO and LUMO of cis butadiene and determine its symmetry (symmetric or anti-symmetric) with respect to plane.

''There are two ways to do this in GaussView. One is: Select '''Edit→MOs'''. Select the HOMO and the LUMO from the MO list (highlights it yellow). Click the button '''Visualise''' (not Calculation), then '''Update'''. Alternately, having calculated the surface for this orbital, you can display it in the main GaussView window for the molecule, from the '''Results→Surfaces''' menu. Select '''Surface Actions→Show Surface'''. Having displayed the surface this way, you can also select '''View→Display Format→Surface''', and change '''Solid''' to '''Mesh'''.''

 

ii) Computation of the Transition State geometry for the prototype reaction and an examination of the nature of the reaction path.

[[Image:mb_da3.jpg |right|thumb|]]

The transition structure has an envelope type structure, which maximizes the overlap between the ethylene π orbitals and the π system of butadiene. One way to obtain the starting geometry is to build the bicyclo system (b) and then remove the -CH2-CH2- fragment. One must then guess the interfragment distance (dashed lines) and optimize the structure, but use any method you wish, based on the tutorial above, to characterise the transition structure. Confirm you have obtained a transition structure for the Diels Alder reaction!

[[Image:mb_da4.jpg |right|thumb|guessing the transition structure]]

Once you have obtained the correct structure, plot the HOMO as in (i). Rotate the molecule so that the symmetry and nodal properties of the system can be interpreted, and save a copy of the image.

 

(iii) To Study the regioselectivity of the Diels Alder Reaction

Cyclohexa-1,3-diene '''1''' undergoes facile reaction with maleic anhydride '''2''' to give primarily the endo adduct. The reaction is supposed to be kinetically controlled so that the exo transition state should be higher in energy.

[[Image:Bearpark_pic_edit_by_jm906.JPG |right|thumb|regioslectivity]]

Locate the transition structures for both 3 and 4. Compare the energies of the endo and exo forms.

Measure the bond lengths of the partly formed σ C-C bonds and the other C-C distances. Make a sketch with the important bond lengths. Measure the orientation, (C-C through space distances between the -(C=O)-O-(C=O)- fragment of the maleic anhydride and the C atoms of the “opposite” -CH2-CH2- for the exo and the “opposite” -CH=CH- for the endo). The structure must be a compromise between steric repulsions of the -CH2-CH2- fragment and the maleic anhydride for the exo versus secondary orbital interactions between the π systems of -CH=CH- and -(C=O)-O-(C=O)- fragment for the endo.

Plot the HOMO as in the previous exercise. Examine carefully the nodal properties of the HOMO between the -(C=O)-O-(C=O)- fragment and the remainder of the system. What can you conclude about the so called “secondary orbital overlap effect”?

=== Suggested Discussion ===

Use this template as a guide. Screen images can be saved from the GaussView '''File''' menu.

''For cis butadiene'': 
Plot the HOMO and LUMO and determine the symmetry (symmetric or anti-symmetric) with respect to the plane.

''For the ethylene+cis butadiene transition structure'': 
Sketch HOMO and LUMO, labeling each as symmetric or anti symmetric.

Show the geometry of the transition structure, including the bond-lengths of the partly formed σ C-C bonds.

What are typical sp3 and sp2 C-C bondlengths? What is the van der Waals radius of the C atom? What can you conclude about the C-C bondlength of the partly formed σ C-C bonds in the TS.

Illustrate the vibration that corresponds to the reaction path at the transition state.
Is the formation of the two bonds synchronous or asynchronous?
How does this compare with the lowest positive frequency?

Is the HOMO at the transition structure '''s''' or '''a'''?

Which MOs of butadiene and ethylene have been used to form this MO?
Explain why the reaction is allowed.

''For the cyclohexa-1,3-diene reaction with maleic anhydride'': 
Give the relative energies of the exo and endo transition structures.
Comment on the structural difference between the endo and exo form. Why do you think that the exo form could be more strained?
Examine carefully the nodal properties of the HOMO between the -(C=O)-O-(C=O)- fragment and the remainder of the system. What can you conclude about the so called “secondary orbital overlap effect”?
(There is some discussion of this in Ian Fleming's book 'Frontier Orbitals and Organic Chemical Reactions').

''Further discussion'': 
What effects have been neglected in these calculations of Diels Alder transition states?

Look at published examples and investigate further if you have time.
(e.g. {{DOI|10.1021/jo0348827}})

----
See also: [[Mod:timetable|Timetable]],[[Mod:lectures|Intro lecture]], [[mod:programs|Programs]], [[mod:organic|Module 1]], [[Mod:inorganic|Module 2]], [[Mod:phys3|Module 3]]

© 2008-2011, Imperial College London

Mod:phys3

2012-10-29T18:02:23Z

Lmt09: /* Optimizing the "Chair" and "Boat" Transition Structures */

See also:[[Mod:timetable|Timetable]], [[Mod:lectures|Intro lecture]], [[mod:programs|Programs]], [[mod:organic|Module 1]], [[Mod:inorganic|Module 2]], [[Mod:phys3|Module 3]], [http://www.gaussian.com/g_tech/gv5ref/gv5ref_toc.htm Gaussian Online User Manual] | [http://faculty.ycp.edu/~jforesma/educ/visual/index.html Visualization Tutorials]
= Module 3 =

In this set of computational experiments, you will characterise transition structures on potential energy surfaces for the Cope rearrangement and Diels Alder cycloaddition reactions.

There are two parts:
a) tutorial material: how to use the programs and methods,
b) more challenging examples, with guidelines but fewer explicit instructions.



In the second year physical chemistry laboratory, you may have carried out dynamics calculations using model potential energy surfaces to explore transition states. In that computational experiment, the total energy could quickly be calculated for different geometries of a triatomic system using an analytical function of the atomic coordinates (for more information, see for example [http://books.google.com/books?id=T8IZ1aa_FRkC&pg=RA1-PA36&lpg=RA1-PA36&dq=%22lake+eyring%22&source=web&ots=OXY00lSZ7D&sig=Ld_MTNwNjUDNGzB_5w1IxaMBMPU&hl=en&sa=X&oi=book_result&resnum=7&ct=result here] and [http://www.rsc.org/ejarchive/DC/1979/DC9796700007.pdf here]).

In this experiment, you will be studying transition structures in larger molecules. There are no longer fitted formulae for the energy, and the molecular mechanics / force field methods that work well for structure determination cannot be used (in general) as they do not describe bonds being made and broken, and changes in bonding type / electron distribution. (This is the main difference from Module 1). Instead, we use molecular orbital-based methods, numerically solving the Schrodinger equation, and locating transition structures based on the local shape of a potential energy surface. As well as showing what transition structures look like, reaction paths and barrier heights can also be calculated.

==The Cope Rearrangement Tutorial==



''This part of the module is described as a 'tutorial' because it's an introduction to various computational techniques for locating transition structures on potential energy surfaces. It's different to the GaussView tutorial you may have worked through earlier: it's an exercise where you're given specific instructions, see if you can follow them, and also whether there are problems or better ways of carrying the exercise out. Please include this part in your write-up. Marks will be given for correct answers, the documentation showing how you got these, discussion, and how you went about solving any problems you encountered.''

In this tutorial we will use the Cope rearrangement of 1,5-hexadiene as an example of how to study a chemical reactivity problem.

Your objectives are to locate the low-energy minima and transition structures on the C6H10 potential energy surface, to determine the preferred reaction mechanism.

[[Image:pic1.jpg|right|thumb|Cope rearrangement]]

This [3,3]-sigmatropic shift rearrangement has been the subject of numerous experimental and computational studies (e.g. Houk et al. {{DOI|10.1021/ja00101a078}}), and for a long time its mechanism (concerted, stepwise or dissociative) was the subject of some controversy. Nowadays it is generally accepted that the reaction occurs in a concerted fashion via either a "chair" or a "boat" transition structure, with the "boat" transition structure lying several kcal/mol higher in energy. The B3LYP/6-31G* level of theory has been shown to give activation energies and enthalpies in remarkably good agreement with experiment. In this tutorial we will show how these can be calculated using Gaussian.

 

{| align="center"
|- align="center"
|
[[Image:pic2a.jpg]]
|
[[Image:pic2b.jpg]]
|- align="center"
| align="center" | ''Chair Transition State''
| align="center" | ''Boat Transition State''
|}

 

===Optimizing the Reactants and Products===

'' In this section you will learn how to optimize a structure, symmetrize it to find its point group, calculate and visualize vibrational frequencies and correct potential energies in order to compare them with experimental values. It is assumed that you are already familiar with using the builder in GaussView. ''

 

(a) Using GaussView, draw a molecule of 1,5-hexadiene with an "anti" linkage (aproximately a.p.p conformation) for the central four C atoms . Clean the structure using the '''Clean''' function under the '''Edit''' menu.

Now we will optimize the structure at the '''HF/3-21G''' level of theory. Select '''Gaussian''' under the '''Calculate''' menu, click on the '''Job Type''' tab and choose '''Optimization'''. The default method should already be Hartree Fock and the default basis set is 3-21G, so there should be no need to change these. You can check this by clicking on the '''Method''' tab. Change the %mem under the '''Link 0''' tab to 250 MB (though this could be increased to 500 MB). Submit the job by clicking on the '''Submit''' button at the bottom of the window and give the job a meaningful name (e.g. react_anti).

When the job has finished, you will be asked if you want to open a file. Select '''Yes''' and choose the checkpoint (chk) file with the name of the job you have just run (e.g. react_anti.chk). This checkpoint file is a binary file that stores data calculated by Gaussian. The name of the chk file should have been assigned by default, but by default, this file will be created in the C:\Windows\G03\Scratch folder. Once the file has been opened, click on the '''Summary''' button under the '''Results''' menu and make a note of the energy.


Does your final structure have symmetry? Select '''Symmetrize''' under the '''Edit''' menu (note that sometimes it is necessary to relax the search criteria under the '''Point Group''' menu). Make a note of the point group.

 

(b) Now draw another molecule of 1,5-hexadiene with a "gauche" linkage for the central four C atoms. Would you expect this structure to have a lower or a higher energy than the anti structure you have just optimized? Optimize the structure at the '''HF/3-21G''' level of theory and compare your final energy with that obtained in (a). Again, check if the molecule has symmetry and make a note of the point group.

 

(c) Normally, calculated activation energies and enthalpies use the lowest energy conformation of a reactant molecule as a reference. Based on your results from above, try to predict what the lowest energy conformation of 1,5-hexadiene might be. Test out your hypothesis by drawing the structure and optimizing it.

 

(d) A table containing the low energy conformers of 1,5-hexadiene and their point groups is shown in [[Mod:phys3#Appendix 1|Appendix 1]]. Compare the structures that you have optimized with those in the table and see if you can identify your structure.

 

(e) Draw the Ci ''anti2'' conformation of 1,5-hexadiene (unless you have already located it). Optimize it at the '''HF/3-21G''' level of theory and make sure it has Ci symmetry. Compare your final energy to the one given in the table.


 

(f) When you are happy that your structure is the same as the one in the table, reoptimize it at the '''B3LYP/6-31G*''' level (6-31G* is equivalent to 6-31G(d) by selecting '''DFT''' under the '''Method''' menu and '''B3LYP''' from the box with the functionals on the right-hand side. Now select '''Link 0''' and change the name of the chk file to the name of the DFT optimization that you are about to run. Note that it is always advisable to do this when re-using or modifying existing structures to ensure that the original chk file is not overwritten. Run the job and make a note of the energy. Now compare the final structures from the '''HF/3-21G''' calculation with that at the higher level of theory. How much does the overall geometry change?

 

(g) The final energies given in the output file represent the energy of the molecule on the bare potential energy surface. To be able to compare these energies with experimentally measured quantities, they need to include some additional terms, which requires a frequency calculation to be carried out. The frequency calculation can also be used to characterize the critical point, i.e. to confirm that it is a minimum in this case: that all vibrational frequencies are real and positive.

Starting from your optimized B3LYP/6-31G* structure, run a frequency calculation at the same level of theory. You can do this by selecting '''Frequency''' under the '''Job Type''' tab. Ensure that the method is still correctly specified under the '''Method''' tab (''caution: on Windows, sometimes 'scrf=(solvent=water,check)' is incorrectly added!'') and then change the name of the chk file under the '''Link 0''' tab to the name of the frequency job that you are about to run. Run the job. Once the job has finished, open the log file this time. Select '''Vibrations''' under the '''Results''' menu. A list of all the vibrational frequencies modes should appear. Check that there are no imaginary frequencies, only real ones. You can visualize some of these vibrations under this menu and simulate the infrared spectrum.


Now, select '''View File''' under the '''Results''' menu and open the output file in the visualizer. Scroll down to the section beginning '''Thermochemistry'''. Under the vibrational temperatures a list of energies should be printed. Make a note of (i) the sum of electronic and zero-point energies, (ii) the sum of electronic and thermal energies, (iii) the sum of electronic and thermal enthalpies, and (iv) the sum of electronic and thermal free energies. The first of these is the potential energy at 0 K including the zero-point vibrational energy (E = Eelec + ZPE), the second is the energy at 298.15 K and 1 atm of pressure which includes contributions from the translational, rotational, and vibrational energy modes at this temperature (E = E + Evib + Erot + Etrans), the third contains an additional correction for RT (H = E + RT) which is particularly important when looking at dissociation reactions, and the last includes the entropic contribution to the free energy (G = H - TS). It is important to make sure that you select the correct energy/enthalpy term to compare to your experimental values. Note that these corrections can also be calculated at other temperatures using the '''Freq=ReadIsotopes''' option in Gaussian, If you have time, try re-calculate these quantities at 0 K as shown in the [[mod:gv_advanced | Advanced GaussView Tutorial]].

 

===Optimizing the "Chair" and "Boat" Transition Structures ===

'' In this section you will learn how to set up a transition structure optimization (i) by computing the force constants at the beginning of the calculation, (ii) using the redundant coordinate editor, and (iii) using QST2. You will also visualize the reaction coordinate and run the IRC (Intrinisic Reaction Coordinate) and calculate the activation energies for the Cope rearrangement via the "chair" and "boat" transition structures. ''

 

The "chair" and "boat" transition structures for the Cope rearrangement are shown in [[Mod:phys3#Appendix 2|Appendix 2]]. Both consist of two C3H5 allyl fragments positioned approximately 2.2 Å apart, one with C2h symmetry and the other with C2v symmetry.

 

(a) Draw an allyl fragment (CH2CHCH2) and optimize it using the '''HF/3-21G''' level of theory. Your structure should look like one half of the transition structures shown below.

Now open a new GaussView window by going to the '''File''' menu and selecting '''New''' and then '''Create MolGroup'''. Copy the optimized allyl structure from the first calculation by selecting '''Copy''' under the '''Edit''' menu, and then paste it twice into the new window by selecting '''Paste''' and then '''Append Molecule'''. Now orient the two fragments so that they look roughly like the chair transition state below by using the '''Shift Alt keys + Left Mouse button''' to translate one fragment with respect to the other and the '''Alt key + Left Mouse button''' to rotate it. The distance between the terminal ends of the allyl fragments should be approximately 2.2 Å apart. Save this structure to a Gaussian input file with a meaningful name (e.g. chair_ts_guess).

We are now going to optimize this transition state manually in two different ways. Transition state optimizations are more difficult than minimizations because the calculation needs to know where the negative direction of curvature (i.e. the reaction coordinate) is. If you have a reasonable guess for your transition structure geometry, then normally the easiest way to produce this information is to compute the force constant matrix (also known as the Hessian) in the first step of the optimization which will then be updated as the optimization proceeds. This is what we will try to do in the next section. However, if the guess structure for the transition structure is far from the exact structure, then this approach may not work as the curvature of the surface may be significantly different at points far removed from the transition structure. In some cases, a better transition structure can be generated by freezing the reaction coordinate (using '''Opt=ModRedundant''' and minimizing the rest of the molecule. Once the molecule is fully relaxed, the reaction coordinate can then be unfrozen and the transition state optimization is started again. One advantage of doing this, is that it may not be necessary to compute the whole Hessian once this has been done, and just differentiating along the reaction coordinate might give a good enough guess for the initial force constant matrix. This can save a considerable amount of time in cases where the force constant calculation is expensive.

 

(b) Use Hartree Fock and the default basis set 3-21G for parts (b) to (f).

Create a new MolGroup ('''File → New → Create MolGroup''') and copy and paste your guess structure into the window. Now set up a Gaussian optimization for a transition state. Go to the '''Gaussian''' menu under '''Calculate''' and click on the '''Job Type''' tab. Select '''Opt+Freq''' and then change '''Optimization to a Minimum''' to '''Optimization to a TS (Berny)'''. Choose to calculate the force constants '''Once''' and in the Additional keyword box at the bottom, type '''Opt=NoEigen'''. The latter stops the calculation crashing if more than one imaginary frequency is detected during the optimization which can often happen if the guess transition structure is not good enough. Submit the job. If the job completes successfully, you should have optimized to the structure shown in [[Mod:phys3#Appendix 2|Appendix 2]] and the frequency calculation should give an imaginary frequency of magnitude 818 cm-1. Animate the vibration and ensure that it is the one corresponding to the Cope rearrangement.

 

(c) Now we will try optimizing the transition structure again using the frozen coordinate method. Create a new MolGroup ('''File → New → Create MolGroup''') and copy and paste your guess structure into the window again. Now select '''Redundant Coord Editor''' from the '''Edit''' menu. Click on the highlighted file icon at the top left-hand corner (Create a New Coordinate) and a line should appear below saying '''Add Unidentified (?, ?, ?, ?)'''. Now go back to the GaussView window and select two of the terminal carbons from the allyl fragments which form/break a bond during the rearrangement. Return to the coordinate editor and select '''Bond''' instead of '''Unidentified''' and select '''Freeze Coordinate''' instead of '''Add'''. Now click on the icon again to generate another coordinate. This time select the opposite two terminal atoms and again select '''Bond''' and '''Freeze Coordinate'''. Click OK. Now set up the optimization as if it were a minimum and you should see the option '''Opt=ModRedundant''' already included in the input line. Submit the job.

'''Note:''' GaussView allows you to produce an input file with the frozen coordinate specified as e.g. <tt>B 5 1 2.200000 F</tt>. Unfortunately, a recent update to the Gaussian program means it does not recognise this syntax, and just ignores this line. This means that the coordinate ends up being optimised rather than frozen. Therefore do not use this method, but ensure the guess structure has suitable guess transition bond distances(~2.2 Å) using the ''Modify Bond'' tool in GaussView --[[User:Rzepa|Rzepa]] 14:39, 29 October 2012 (UTC)

 

(d) When the job has finished, open the chk file. You should find that the optimized structure looks a lot like the transition you optimized in section (b), except the bond forming/breaking distances are fixed to 2.2 Å. Now we are going to optimize them too. Open the '''Redundant Coord Editor''' from the '''Edit''' menu again and create a new coordinate as before by clicking on the icon, Select one of the bonds that was previously frozen and this time choose '''Bond''' instead of '''Unidentified''' and '''Derivative''' instead of '''Add'''. Repeat the procedure for the other bond. This time you need to set up a transition state optimization but we are not going to calculate the force constants as we did in section (b) (so we leave this option as '''Never'''), instead we will use a normal guess Hessian modified to include the information about the two coordinates we are differentiating along. Change the name of the chk file in '''Link 0''' if you do not want to write over the previous calculation and submit the job. When the calculation has finished, open the chk file, check the bond forming/bond breaking bond lengths and compare the structure to the one you optimized in section (b).

 

(e) Now we will optimize the boat transition structure. We will do this using the '''QST2''' method. In this method, you can specify the reactants and products for a reaction and the calculation will interpolate between the two structures to try to find the transition state between them. You must make sure that your reactants and products are numbered in the same way. Therefore, although our reactants and products are both 1,5-hexadiene, we will need to manually change the numbering for the product molecule so that it corresponds to the numbering obtained if our reactant had rearranged.

e.g.

<center>[[Image:pic3.jpg|200px]]</center>

Open the chk file corresponding to the optimized Ci reactant molecule (''anti2'' in [[Mod:phys3#Appendix 1|Appendix 1]]). Now open a second window and create a new MolGroup. Copy the optimized reactant molecule into the new window. In the same window, now select '''File → New → Add to MolGroup'''. The original molecule should disappear and a green circle should appear at the top left-hand corner with a '''2''' next to it. Clicking on the down arrow by the '''2''' will take you back to the original window and you will see your molecule again. This is how we read multiple geometries into GaussView. Go back to window '''2''', and copy and paste the reactant molecule a second time. This is going to be the product molecule and will be the molecule on which we need to change the numbering. If you now click on the icon showing two molecules side by side, then you can view both molecules simultaneously.

Now go to the '''View''' menu and select '''Labels''' so that you can see the numbering on both structures. Orient the two structures separately so they look something like the following:

{| align="center"
|- align="center"
|
[[Image:pic4a.jpg|200px]]
|
[[Image:pic4b.jpg|200px]]
|- align="center"
| align="center" | ''Reactant''
| align="center" | ''Product''
|}

 

Now click on the product structure. Go to the '''Edit''' menu and select '''Atom List'''. Starting from Atom 1 on the reactant, go through and renumber all the atoms on the Product so that they match the reactant molecule, e.g. for the numbering above you would start by changing atom '''6''' on the product molecule to atom '''1'''. The other atom numbers will update as you do this so make sure you do it in the correct order. At the end, the numbering on your two molecules should correspond to each other in the following way:

 

{| align="center"
|- align="center"
|
[[Image:pic5a.jpg|200px]]
|
[[Image:pic5b.jpg|200px]]
|- align="center"
| align="center" | ''Reactant''
| align="center" | ''Product''
|}

 

Now we will set up the first '''QST2''' calculation. Go to the '''Gaussian''' menu and select '''Job Type''' as '''Opt+Freq''', and optimize to a transition state. This time you will have two options - '''TS (Berny)''' which we used in the previous calculations and '''TS (QST2)'''. Select '''TS (QST2)'''. Submit the job.

You will find that the job fails. To see why, open the chk file you created and view the structure. You will see that it looks a bit like the chair transition structure but more dissociated. In fact when the calculation linearly interpolated between the two structures, it simply translated the top '''allyl''' fragment and did not even consider the possibility of a rotation around the central bonds. It is clear that the QST2 method is never going to locate the boat transition structure if we start from these reactant and product structures.

Now go back to the original input file where you set up your QST2 calculation. We will now modify the reactant and product geometries so that they are closer to the boat transition structure. Click on the reactant molecule first and select the central '''C-C-C-C''' dihedral angle (i.e. '''C2-C3-C4-C5''' for the molecule above) and change the angle to 0o. Then select the inside '''C-C-C''' (i.e. '''C2-C3-C4''' and '''C3-C4-C5''' for the molecule above) and reduce them to 100o. Do the same for the product molecule. Your reactant and product molecules should now look like the following:

 

{| align="center"
|- align="center"
|
[[Image:pic6a.jpg|200px]]
|
[[Image:pic6b.jpg|200px]]
|- align="center"
| align="center" | ''Reactant''
| align="center" | ''Product''
|}

 

Set up the QST2 calculation again, renaming both the chk file under '''Link 0''' and the input file. Run the job again. This time it should converge to the boat transition structure. Check that there is only one imaginary frequency and visualize its motion.

The object of this exercise is to illustrate that although the QST2 method is has some advantages because it is fully automated, it can often fail if your reactants and products are not close to the transition structure. There is another method, the '''QST3''' method, that allows you to input the geometry of a guess transition structure also and this can often be more reliable. If you have time, you can try generating a guess boat transition structure and see if you can get the calculation to converge using the original reactant and product molecules. Remember to check the atom numbers in the transition structure are in the right order.

 

(f) Take a look at your optimized chair and boat transition structures. Which conformers of 1,5-hexadiene do you think they connect? You will find that it is almost impossible to predict which conformer the reaction paths from the transitions structures will lead to. However, there is a method implemented in Gaussian which allows you to follow the minimum energy path from a transition structure down to its local minimum on a potential energy surface. This is called the '''Intrinsic Reaction Coordinate''' or '''IRC''' method. This creates a series of points by taking small geometry steps in the direction where the gradient or slope of the energy surface is steepest.

Open the chk file for one of your optimized chair transition structures. Under the '''Gaussian''' menu, select '''IRC''' under the '''Job Type''' tab. You will be presented with a number of options. The first is to decide whether to compute the reaction coordinate in one or both directions. As our reaction coordinate is symmetrical, we will only choose to compute it in the forward direction. Normally you would do both forward and reverse, either in one job or in two separate jobs. You are also given the option to calculate the force constants once, at every step along the IRC or to read them from the chk file. You would use the latter option if you have previously run a frequency calculation. In this case, to avoid confusion with chk files, we will just recompute them at the beginning of the calculation. (The '''IRCMax''' option can also be specified here. This takes a transition structure as its input, and finds the maximum energy along a specified reaction path, taking into account zero-point energy etc., and produces all the quantities needed for a variational transition state theory calculation. We will leave this unchecked for the purposes of this exercise.) The final option to consider is the number of points along the IRC. The default is '''6''' but this is normally never enough. Let's change this to 50 and see how the calculation progresses. Change the name of the chk file under '''Link 0''' and submit the job. The job will take a while so now is a good time to take a coffee break...

When the IRC calculation has finished, open the chk file with all the intermediate geometries and see how the calculation has progressed. You will find that it hasn't reached a minimum geometry yet. This leaves you three options: (i) you can take the last point on the IRC and run a normal minimization; (ii) you can restart the IRC and specify a larger number of points until it reaches a minimum; (iii) you can redo the IRC specifying that you want to compute the force constants at every step. There are advantages and disadvantages to each of these approaches. Approach (i) is the fastest, but if you are not close enough to a local minimum, you may end up in the wrong minimum. Approach (ii) is more reliable but if too many points are needed, then you can also veer off in the wrong direction after a while and end up at the wrong structure. Approach (iii) is the most reliable but also the most expensive and is not always feasible for large systems. You can try any or all of these approaches and see which conformation you end up in.

 

(g) Finally we need to calculate the activation energies for our reaction via both transition structures. To do this we will need to reoptimize the chair and boat transition structures using the '''B3LYP/6-31G*''' level of theory and to carry out frequency calculations. You can start from the HF/3-21G optimized structures. Once the calculations have converged, compare both the geometries and the difference in energies between the reactants and transition states at the two levels of theory. What you should find is that the geometries are reasonably similar, but the energy differences are markedly different.

As a consequence of this, it is often more computational efficient to map the potential energy surface using the low level of theory first and then to reoptimize at the higher level as we have done in this exercise.

The experimental activation energies are 33.5 ± 0.5 kcal/mol via the chair transition structure and 44.7 ± 2.0 kcal/mol via the boat transition structure at 0 K. If you take the values computed at 0 K, how close are they to the experimental values? You can also find the energies with thermal correction at 298.15 K under the Thermochemistry data in the output file. If you have time, you can recompute them at higher temperature. Alternatively, you can use the utility program '''FreqChk''' to obtain energies at a different temperature. This only requires the chk file from a frequency calculation and allows you to retrieve frequency and thermochemistry data as well as calculating them with an alternate temperature, pressure, scale factor, and/or isotope substitutions. The '''FreqChk''' utility program can be accessed from '''Gaussian03W'''. Launch '''Gaussian03W'''. Select '''utilities''' from the menu and click on '''FreqChk''' to launch the utility program. You will be prompted for a chk file. Select your chk file from the C:\G03W\Scratch directory and follow the instructions from this [http://www.gaussian.com/g_tech/g_ur/u_freqchk.htm web link] to proceed.

 

=== Appendix 1 ===

 

{| border="1" cellpadding="5"
|- align="center"
| width="150" | '''Conformer'''
| width="150" | '''Structure'''
| width="100" | '''Point Group'''
| width="200" | '''Energy/Hartrees HF/3-21G'''
| width="200" | '''Relative Energy/kcal/mol'''
|- align="center"
| ''gauche1''
|
[[Image:gauche1.jpg|150px]]
| C2
| -231.68772
| 3.10
|- align="center"
| ''gauche2''
|
[[Image:gauche2.jpg|150px]]
| C2
| -231.69167
| 0.62
|- align="center"
| ''gauche3''
|
[[Image:gauche3.jpg|150px]]
| C1
| -231.69266
| 0.00
|- align="center"
| ''gauche4''
|
[[Image:gauche4.jpg|150px]]
| C2
| -231.69153
| 0.71
|- align="center"
| ''gauche5''
|
[[Image:gauche5.jpg|150px]]
| C1
| -231.68962
| 1.91
|- align="center"
| ''gauche6''
|
[[Image:gauche6.jpg|150px]]
| C1
| -231.68916
| 2.20
|- align="center"
| ''anti1''
|
[[Image:anti1.jpg|150px]]
| C2
| -231.69260
| 0.04
|- align="center"
| ''anti2''
|
[[Image:anti2.jpg|150px]]
| Ci
| -231.69254
| 0.08
|- align="center"
| ''anti3''
|
[[Image:anti3.jpg|150px]]
| C2h
| -231.68907
| 2.25
|- align="center"
| ''anti4''
|
[[Image:anti4.jpg|150px]]
| C1
| -231.69097
| 1.06
|}

 

=== Appendix 2 ===

{| cellpadding="5"
|- align="center"
|
[[Image:appendix2a.jpg|300px]]
|- align="center"
| ''C2h Chair Transition State''
|- align="center"
|- align="center"
|- align="center"
|
[[Image:appendix2b.jpg|300px]]
|- align="center"
| ''C2v Boat Transition State''
|}

 

=== Results Table ===

 

''' Summary of energies (in hartree) '''

{| border="1" cellspacing="1" cellpadding="10"
|-
| ''' '''
!colspan="3"|'''HF/3-21G'''
!colspan="3"|'''B3LYP/6-31G*'''
|-
| ''' '''
| width="125" align="center" | '''Electronic energy'''
| width="125" align="center" | '''Sum of electronic and zero-point energies'''
| width="125" align="center" | '''Sum of electronic and thermal energies'''
| width="125" align="center" | '''Electronic energy'''
| width="125" align="center" | '''Sum of electronic and zero-point energies'''
| width="125" align="center" | '''Sum of electronic and thermal energies'''
|-
| ''' '''
| ''' '''
| width="125" align="center" | '''at 0 K'''
| width="125" align="center" | '''at 298.15 K'''
| ''' '''
| width="125" align="center" | '''at 0 K'''
| width="125" align="center" | '''at 298.15 K'''
|-
| '''Chair TS'''
| align="center" | -231.619322
| align="center" | -231.466705
| align="center" | -231.461346
| align="center" | -234.556983
| align="center" | -234.414919
| align="center" | -234.408998
|-
| '''Boat TS'''
| align="center" | -231.602802
| align="center" | -231.450929
| align="center" | -231.445300
| align="center" | -234.543093
| align="center" | -234.402340
| align="center" | -234.396006
|-
| '''Reactant (''anti2'')'''
| align="center" | -231.692535
| align="center" | -231.539539
| align="center" | -231.532566
| align="center" | -234.611710
| align="center" | -234.469203
| align="center" | -234.461856
|}

 *1 hartree = 627.509 kcal/mol 

''' Summary of activation energies (in kcal/mol) '''

{| border="1" cellspacing="1" cellpadding="10"
|-
| ''' '''
| width="125" align="center" | '''HF/3-21G'''
| width="125" align="center" | '''HF/3-21G'''
| width="125" align="center" | '''B3LYP/6-31G*'''
| width="125" align="center" | '''B3LYP/6-31G*'''
| width="125" align="center" | '''Expt.'''
|-
| ''' '''
| width="125" align="center" | '''at 0 K '''
| width="125" align="center" | '''at 298.15 K'''
| width="125" align="center" | '''at 0 K'''
| width="125" align="center" | '''at 298.15 K'''
| width="125" align="center" | '''at 0 K'''
|-
| '''ΔE (Chair)'''
| align="center" | 45.70
| align="center" | 44.69
| align="center" | 34.06
| align="center" | 33.17
| align="center" | 33.5 ± 0.5
|-
| '''ΔE (Boat)'''
| align="center" | 55.60
| align="center" | 54.76
| align="center" | 41.96
| align="center" | 41.32
| align="center" | 44.7 ± 2.0
|}

 

==The Diels Alder Cycloaddition==

In this exercise, you will characterise transition structures using any of the methods described above in the tutorial: the choice is up to you. In addition, you will look at the shape of some of the molecular orbitals. To help you structure your report, there is a data/discussion sheet at the end of this section.

[[Image:mb_da1.jpg |right|thumb|Diels Alder cycloaddition]]
The Diels Alder reaction belongs to a class of reactions known as pericyclic reactions. The π orbitals of the dieneophile are used to form new σ bonds with the π orbitals of the diene. Whether or not the reactions occur in a concerted stereospecific fashion ('''allowed''') or not ('''forbidden''') depends on the number of π electrons involved. In general the HOMO/LUMO of one fragment interacts with the HOMO/LUMO of the other reactant to form two new bonding and anti-bonding MOs. The nodal properties allow one to make predictions according to the following rule:

''If the HOMO of one reactant can interact with the LUMO of the other reactant then the reaction is '''allowed'''.''

''The HOMO-LUMO can only interact when there is a significant overlap density. If the orbitals have different symmetry properties then no overlap density is possible and the reaction is '''forbidden'''.''

If the dieneophile is substituted, with substituents that have π orbitals that can interact with the new double bond that is being formed in the product, then this interaction can stabilise the regiochemistry (i.e. head to tail versus tail to head) of the reaction. In this exercise you will study the nature of the transition structure for the Diels Alder reaction, both for the prototypical reaction and for the case where both diene and dieneophile carry substituents, and where secondary orbital effects are possible. Clearly, the factors that control the nature of the transition state are quantum mechanical in origin and thus we shall use methods based upon quantum chemistry.

 

Shown on the right is a diagram of the transition state for the Diels-Alder reaction between ethylene and butadiene. The ethylene approaches the cis form of butadiene from above.
[[Image:mb_da2.jpg |right|thumb|Ethylene+Butadiene cycloaddition]]

Before beginning our quantitative study, it is helpful to discuss the interaction of the π orbitals in a simple qualitative way.

''You will confirm some of these considerations in your computations.''

The principal orbital interactions involve the π/ π* orbitals of ethylene and the HOMO/LUMO of butadiene. It is referred to as [4s + 2s] since one has 4 π orbitals in the π system of butadiene. The orbitals of ethylene and butadiene and ethylene can be classified as symmetric '''s''' or anti-symmetric '''a''' with respect to the plane of symmetry shown.

The HOMO of ethylene and the LUMO of butadiene are both '''s''' (symmetric with respect to the reflection plane) and the LUMO of ethylene and the HOMO of butadiene are both '''a'''. Thus it is the HOMO-LUMO pairs of orbital that interact, and energetically, the HOMO of the resulting adduct with two new σ bonds is '''a'''.

 

=== Exercise ===

Use the the AM1 semi-empirical molecular orbital method for these calculations (to start with).

i) Use GaussView to build cis butadiene, and optimize the geometry using Gaussian. Plot the HOMO and LUMO of cis butadiene and determine its symmetry (symmetric or anti-symmetric) with respect to plane.

''There are two ways to do this in GaussView. One is: Select '''Edit→MOs'''. Select the HOMO and the LUMO from the MO list (highlights it yellow). Click the button '''Visualise''' (not Calculation), then '''Update'''. Alternately, having calculated the surface for this orbital, you can display it in the main GaussView window for the molecule, from the '''Results→Surfaces''' menu. Select '''Surface Actions→Show Surface'''. Having displayed the surface this way, you can also select '''View→Display Format→Surface''', and change '''Solid''' to '''Mesh'''.''

 

ii) Computation of the Transition State geometry for the prototype reaction and an examination of the nature of the reaction path.

[[Image:mb_da3.jpg |right|thumb|]]

The transition structure has an envelope type structure, which maximizes the overlap between the ethylene π orbitals and the π system of butadiene. One way to obtain the starting geometry is to build the bicyclo system (b) and then remove the -CH2-CH2- fragment. One must then guess the interfragment distance (dashed lines) and optimize the structure, but use any method you wish, based on the tutorial above, to characterise the transition structure. Confirm you have obtained a transition structure for the Diels Alder reaction!

[[Image:mb_da4.jpg |right|thumb|guessing the transition structure]]

Once you have obtained the correct structure, plot the HOMO as in (i). Rotate the molecule so that the symmetry and nodal properties of the system can be interpreted, and save a copy of the image.

 

(iii) To Study the regioselectivity of the Diels Alder Reaction

Cyclohexa-1,3-diene '''1''' undergoes facile reaction with maleic anhydride '''2''' to give primarily the endo adduct. The reaction is supposed to be kinetically controlled so that the exo transition state should be higher in energy.

[[Image:Bearpark_pic_edit_by_jm906.JPG |right|thumb|regioslectivity]]

Locate the transition structures for both 3 and 4. Compare the energies of the endo and exo forms.

Measure the bond lengths of the partly formed σ C-C bonds and the other C-C distances. Make a sketch with the important bond lengths. Measure the orientation, (C-C through space distances between the -(C=O)-O-(C=O)- fragment of the maleic anhydride and the C atoms of the “opposite” -CH2-CH2- for the exo and the “opposite” -CH=CH- for the endo). The structure must be a compromise between steric repulsions of the -CH2-CH2- fragment and the maleic anhydride for the exo versus secondary orbital interactions between the π systems of -CH=CH- and -(C=O)-O-(C=O)- fragment for the endo.

Plot the HOMO as in the previous exercise. Examine carefully the nodal properties of the HOMO between the -(C=O)-O-(C=O)- fragment and the remainder of the system. What can you conclude about the so called “secondary orbital overlap effect”?

=== Suggested Discussion ===

Use this template as a guide. Screen images can be saved from the GaussView '''File''' menu.

''For cis butadiene'': 
Plot the HOMO and LUMO and determine the symmetry (symmetric or anti-symmetric) with respect to the plane.

''For the ethylene+cis butadiene transition structure'': 
Sketch HOMO and LUMO, labeling each as symmetric or anti symmetric.

Show the geometry of the transition structure, including the bond-lengths of the partly formed σ C-C bonds.

What are typical sp3 and sp2 C-C bondlengths? What is the van der Waals radius of the C atom? What can you conclude about the C-C bondlength of the partly formed σ C-C bonds in the TS.

Illustrate the vibration that corresponds to the reaction path at the transition state.
Is the formation of the two bonds synchronous or asynchronous?
How does this compare with the lowest positive frequency?

Is the HOMO at the transition structure '''s''' or '''a'''?

Which MOs of butadiene and ethylene have been used to form this MO?
Explain why the reaction is allowed.

''For the cyclohexa-1,3-diene reaction with maleic anhydride'': 
Give the relative energies of the exo and endo transition structures.
Comment on the structural difference between the endo and exo form. Why do you think that the exo form could be more strained?
Examine carefully the nodal properties of the HOMO between the -(C=O)-O-(C=O)- fragment and the remainder of the system. What can you conclude about the so called “secondary orbital overlap effect”?
(There is some discussion of this in Ian Fleming's book 'Frontier Orbitals and Organic Chemical Reactions').

''Further discussion'': 
What effects have been neglected in these calculations of Diels Alder transition states?

Look at published examples and investigate further if you have time.
(e.g. {{DOI|10.1021/jo0348827}})

----
See also: [[Mod:timetable|Timetable]],[[Mod:lectures|Intro lecture]], [[mod:programs|Programs]], [[mod:organic|Module 1]], [[Mod:inorganic|Module 2]], [[Mod:phys3|Module 3]]

© 2008-2011, Imperial College London

Guide to Creating ONIOM input files for biomolecules

2012-09-05T11:07:58Z

Lmt09: /* Producing a Gaussian Input File */

==Overview==

This guide provides a step by step process to create ONIOM input files for biomolecules from a structure file in the [http://www.pdb.org/pdb/home/home.do Protein Data Bank]. Emphasis will be placed on the use of [https://wiki.ch.ic.ac.uk/wiki/index.php?title=AMBER AMBER] as the low level method and how to obtain parameters for any non-standard residues.

==Creating Standarized .pdb Files==

The first step is to select a .pdb file from the Protein Data Bank that is of high enough resolution to allow atomistic calculations to be produced. The relevant data for determining this is shown on the right hand side under experimental details. The two values to look at are the Resolution[Å] and R-Value, which both should be as low as possible. Having determined a suitable structure, download the suitable test pdb file (usually found in the download files drop-down menu in the top right corner).

In Gaussview select File→Open and choose options. Change the drop-down box "Add Hydrogens:" to '''No''' and, if you wish to remove water molecules, '''check the box''' "Skip Water Molecules." When the file opens up there may be a number of secondary structures present labelled A, B, C etc. In general we require only one so any extras can be removed using Edit→PDB Secondary Structure... and deleting those which are not required. This was then saved as a .pdb file.

{| class="wikitable" border="1"
|-
!The Secondary Structure Editor
|-
|An example of the secondary structure editor is shown below. If we were interested in obtaining structure A only then all that is required is to highlight chains B to D and Edit→Delete→Delete Selected Secondary Structures. The numbered residues such as Helix and Sheet that do not belong to A are automatically removed so if you remove these separately you may end up removing residues from the structure you wish to keep.
[[Image:SSE_PDB_ONIOM_WIKI.png|800 px|alt=PDB Secondary Structure Editor]]
|}

===Residue Names and Protonation States===
Within the .pdb file the fourth column corresponds to the residue name. This name will be used to define the protonation state of the residue, which is currently specified as a default value. In order to determine the protonation state it is possible to use either [http://www.poissonboltzmann.org/pdb2pqr/ PROPKA] or [http://biophysics.cs.vt.edu/ H++]. Once the protonation states have been determined the residue names can be changed to reflect this. (Note: Parameters for non-standard residues calculated later may be included to improve the accuracy).

===Chromophore Structure===
Now a standardized .pdb file of the whole protein has been created the next step is to obtain a .pdb file of the non-standard residue. To do this open the .pdb file we have just saved using a text editor such as vim and remove all lines that are not atoms from the region we intend to include in this residue. It is important here to consider exactly what this consists of here as any problems at this stage are normally not highlighted until much later in the process and will require returning to this point. The region specified here is not the same as that of the ONIOM model region or even the protein chromophore, it is simply so that non-standard residues are defined in the AMBER program. The two important points are that this region must:
# Include the non-standard residue that requires parameterization.
# Is connected to the rest of the protein through standard N or C amino terminations.

The second point may require some elaboration. Some non-standard residues are a modified standard residue, such as that in PYP which is a cystine residue with ''p''-coumaric acid group on the sulphur instead of a thiol. It is tempting to specify the chromophore as just the ''p''-coumaric acid group, however, this causes problems later in defining the parameters for the cystine residue and so the cystine group must also be included in the chromophore region. This joins to the rest of the protein through standard amino acid N and C bonds and so this is all that is needs to be included.

This structure is then saved as a .pdb file and opened in Gaussview. Hydrogens were then added to the residue except where the residue will join to the protein structure. Again be sure of the protonation at this stage as any mistakes will require returning to this point. Check particularly the multiplicity is correct. Save this as a .pdb file and inspect it to ensure that the newly added hydrogens have the same pdb residue name and number as the other atoms, and that their atom numbers follow on and are consistent with connectivity. Also remove any extra TER lines other than the one at the bottom (if there is one). To ensure that this is absolutely correct it may be worth opening this in Gaussview and re-saving it, making sure the correct connectivity is shown.

==Obtaining AMBER Library File of the Chromophore==

We now have two .pdb files, one of the whole protein and one of the non-standard residue region. The next step is to create an AMBER library file of this non-standard residue. Leap, an AMBERTools program, will be used and this requires us to determine three pieces of information for the non-standard residue:
#Connectivity
#AMBER atom types
#Partial charges

Leap can be opened using the command:
<pre>
xleap -s -f /apps/ambertools/amber11/dat/leap/cmd/leaprc.ff03 &
</pre>
This command will use the leap.ff03 set of parameters although any other AMBER parameters could be used depending on the system under study. If this doesn't do anything you probably need to load ambertools:
<pre>
moduleload ambertools
</pre>
Having opened Leap the non-standard residue .pdb file can be loaded using:
<pre>
variable = loadpdb filename
</pre>
where variable is any name you choose and the full pathname must be specified in the filename. Now type:
<pre>
edit variable
</pre>
This brings up a gui where the residue can be visualized. Ensure all atoms are selected and go to Edit→Edit Selected Atoms. This provides a table to be filled with the information specified above. The way to obtain these values will now be explained. A quick sidenote, do not close any Leap x-windows, other than using File→Quit as this will cause the program to crash and any unsaved information to be lost.

====Connectivity====
This is simply achieved by selecting the '''draw''' checkbox in the Leap GUI tool and drawing bonds between the atom centres as desired.

[[Image:Conn_Leap.png|400 px|alt=Drawing bonds with the Leap GUI]]

====AMBER Atom Types====
In order to obtain these open the non-standard residue .pdb file with Gaussview and add methyl groups to the atoms which were previously left with free valences. Save this structure as a .pdb file as we will need it later, however, at this point we only need to go to Edit→Atom List and look at AMBER Type. Copy these across to the Leap table using the PDB Atom Name column to match up Atoms.

====Partial Charges====
This is the most complicated process and requires the use of [[http://q4md-forcefieldtools.org/RED/ R.E.D.-III.4 tools]]. This first uses the modified .pdb file with added methyls to obtain a Gaussian input file using the following command:
<pre>
perl $DIR1/Ante_Red.pl $DIR2/modified_non_standard_residue_file.pdb >> $DIR3/output.txt
</pre>
where $DIR is the relevant pathname. The resulting Gaussian input file can then be run (remember to change memory requirements and checkpoint file locations before submitting). After this has completed the frequency portion was deleted from the log file (this could be removed from the input but is useful for ensuring a minima is obtained) and the log file was copied to Mol_red1.log file in the RED-III directory, ensuring that the filename remains Mol_red1.log. Another file that was output from the above command was a .p2n file. This must be copied to Mol_red1.p2n in the same directory as before, also maintaining Mol_red1.p2n as the filename.

Moving to the RED-III directory now, open Mol_red1.p2n with a text editor and add the following line to exclude the methyl groups from the Partial Charge calculation:
<pre>
REMARK INTRA-MCC 0.0 | 29 30 31 32 33 34 35 36 | R
</pre>
where the numbers correspond to the numbers of the atoms in the methyl groups of the modified non-standard residue. Note that there are two spaces between all the numbers. Below is an example of where it has been placed:

[[Image:p2n_part_char.png|500 px|alt=.p2n file for the calculation of partial charges]]

If necessary change the charge and multiplicity here. Having done this open RED-vIII.4.pl and go to line 4196. Change the variable $DIR to whatever you wish, this is where the output files will be saved to. Create the following jobscript file and run it, although change the directories on line 19 to something useful for you.

<pre>
##################################################################
# REDTOOLS JOBSCRIPT #
# CREATED 08/07/10 #
# LAST MODIFIED 08/07/10 #
# LEE THOMPSON #
##################################################################

#PBS -l ncpus=1
#PBS -l mem=1000mb
#PBS -l walltime=04:00:00
#PBS -joe

module load ambertools
module load gaussian
export GAUSS_SCRDIR=$TMPDIR
echo $GAUSS_SCRDIR
cd $(echo $PBS_O_WORKDIR)
pwd
perl /home/lmt09/SOFTWARE/RED-III.4-Tools-Files/RED-vIII.4.pl > /home/lmt09/PHD_Y2/PYP/1NWZ/PROTONATED/ONIOM/RED_out.log

</pre>

Going into the new directory open Mol-m1-o1-sm.mol2 which contains the partial charges that we seek in the final column. To copy these to the Leap table requires a bit of detective work to match up the atoms. This can be done by opening up the Gaussian log file Mol_red1.log in Gaussview which is labelled in the same order as the .mol2 file with the partial charges. The Gaussview atoms and the atoms in the Leap GUI can then be matched by their positions. These atoms can then be matched to the Leap table by displaying atom names on the Leap GUI using Display→Names. This is also a good time to check consistency of atom types again as if they are different it will cause problems identifying parameters later on. It is also worth checking that the charges sum to an integer value and that '''you have typed them in correctly'''.

----

Having filled in the Leap table go to File→Save and Quit, and then exit the GUI using File→Close. Back at the command line prompt, save the library file using:

<pre>
saveoff variable filename
</pre>
where the variable is the same as used before and filename includes the full pathname. Now exit Leap and go to the .lib file that we have just created. In order for this to be recognised the filename must be uppercase and three or four letters long (although I have not tried to see otherwise). In order to achieve this move it from variable.lib to VAR.lib, where, VAR is a capitalized three letter word of your choice. Now open the file in vi and type:
<pre>
:%s/variable/VAR/g
</pre>
which changes all instances of variable to VAR. We have now created our library file for the non-standard residue.

{| class="wikitable" border="1"
|-
!Justification of Partial Charge Model
|-
|The determination of partial charges is important for the successful use of force field methods, yet the concept of a partial charge is somewhat ambiguous, with several different methods for their determination (see Cramer, C.J., ''Essentials of Computational Chemistry, p309'' for an introduction. The partial charges we use are computed using the restrained ESP method (Cornell ''et al'', ''Journal of the American Chemical Society'', '''1995''', 117, 19, 5179-5197). This is an extension of the ESP method which determines partial charge ''q'' on atom ''k'' by minimizing the difference between:
<math>V_{ESP}(\mathbf{r}) = \sum_{k}^{nuclei}\frac{q_{k}}{\vert \mathbf{r}-\mathbf{r}_{k}\vert}</math><br\>
and the Molecular Electrostatic Potential (MEP):<br\>
<math>V_{MEP}(\mathbf{r}) = \sum_{k}^{nuclei}\frac{Z_{k}}{\vert \mathbf{r}-\mathbf{r}_{k}\vert}-\int \! \Psi (\mathbf{r'}) \frac{1}{\vert \mathbf{r}-\mathbf{r'}\vert} \Psi (\mathbf{r'}) \, \mathrm{d}\mathbf{r'}</math><br\>
for all positions '''r'''. This is computed from a number of points spaced evenly around the Connolly surface of the molecule. ESP is dependent on conformation, however, causing hydrogens in a methyl group for example to have different partial charges. As these are all freely rotating in practice the same partial charges may used for each hydrogen and this is the extension that RESP applies to the ESP method (Bayly ''et al'', ''The Journal of Physical Chemistry'', '''1993''', 97, 40, 10269-10280).
|-
|The main purpose for using this is that AMBER uses RESP for its parm96 (Cornell) parameter set which is the same as that used by Gaussian (derived from HF/6-31G*). Reasons for its use in this force field are that it has been shown to be useful for modeling inter-molecular interactions at short to long range, is convergent with respect to the size of basis set used, resolves to an extent the problems of atoms which do not contribute the Connolly surface and so are ill-defined by the method, as well as having the original advantages of ESP over methods such as Mulliken and Löwdin charges.
|-
|RED (RESP and ESP charge Derive) tools is a series of perl scripts which generate a Gaussian input file which can be run and from which the partial charges derived (Dupradeau ''et al'', ''Physical chemistry chemical physics: PCCP'', '''2010''', 12, 28, 7821-39).
|}

==Producing a Gaussian Input File==

'''NOTE: It is now recommended to add hydrogens using PROPKA or H++ prior to loading into LEAP'''

Having constructed the library file of the non-standard residue we must now construct a .com or .gjf file to run in Gaussian. Initially this will simply be an AMBER calculation, the output of which will be used to determine if we have all the correct parameters and as a starting geometry for the ONIOM calculations. The first step is to reopen Leap using the same command as before. now load in the AMBER library file for the non standard residue using the command:
<pre>
loadoff $DIR/VAR.lib
</pre>
where again $DIR represents the pathname of the file. Typing ''list'' in Leap will display all the library files that have been loaded of which VAR should be one of them. The next stage is to load the .pdb file of the protein that we obtained from Gaussview previously using the command for loading .pdb files shown previously. This should add hydrogens to the structure in accordance to the library files and perhaps a terminal oxygen although never any other heavy atom (this is displayed at the command line). Opening the Leap GUI of the whole protein should reveal the non-standard residue highlighted in the full protein environment. The connection between the protein and the residue must now be determined. There are a number of ways to do this including at the command line but the most successful method so far is to simply draw the bonds in the Leap GUI.

[[Image:Pro_struct.png|500 px|alt=Protein structure in Leap GUI]]

Having drawn the connectivity, go to Unit→Calculate Net Charge to obtain the charge of the protein, which should be an integer. Close the Leap GUI and save the .pdb using the following command:
<pre>
savepdb variable $DIR/filename.pdb
</pre>
We will need the information we have loaded in Leap later so do not close it but for now, look at the .pdb file in the text editor and ensure that there is a terminal oxygen labelled <OXT> at the bottom of the file. If there is not insert ''ATOM 3608 OXT HIE 228 -1.012 21.725 100.791 1.00 0.00'' in the correct place, although the cartesian coordinates, PDB residue name and number and atom number will be different from this example. Also bear in mind that if there are any waters below this then their atom numbers will need changing (use ''grep "WAT" filename1.pdb | awk '{ X=$2; Y=X+1; print "s/"X," "$3" /"Y," "$3" /" }' > sedscript | sed -f sedscript < filename1.pdb > filename2.pdb'' for this). Also check in the .pdb file that the atoms around the chromophore all are of the same residue and do not have differing residue numbers. If this is not the case then it means that the non-standard residue .pdb file was mixed up and you must return to that stage.

Now open the .pdb file in Gaussview and go to Edit→Atom List. Scan through this to ensure that all MM partial charges are present for all atoms other than those in the chromophore residue. If there are any that are not it is because the residue connectivity is wrong so use the bond specification tool to correct this in Gaussview and you should see the MM charges appear as soon as you correct the problem (Hint: The problem atoms will have undefined AMBER atom types (shown as ?) so look at connectivity around these atoms). The MM partial charges can be copied directly from Mol-m1-o1-sm.mol2 into the Gaussview atom list now, although I prefer a second option which I shall explain when I come to it. Now save this as a .com/.gjf file although, because of a bug in Gaussview which causes patial charges to be missing from the input file, you '''must''' save this using Calculate→Gaussian Calculation Setup, chose an AMBER calculation and insert the charge determined earlier and the multiplicity and submit. Select yes when prompted to save the file and then cancel the file execution. You should now have a Gaussian input file in your directory. If you have not inserted the MM partial charges previously copy them from Mol-m1-o1-sm.mol2 and paste them between the AMBER atom type and the PDB information in the Gaussian input file. All that we now require for a complete AMBER calculation is the AMBER parameters for the non-standard residue.

==Getting Non-Standard AMBER Parameters==
If we were to run the Gaussian input file as produced above we would get an error message indicating missing AMBER parameters. Gaussian uses parm96 by default and if any stretches, bends or torsions are present in the non-standard residue but not in the forcefield, then an error message is obtained. In order to obtain them we first need to know what parameters are missing. This can be achieved using:
<pre>
saveamberparm variable xxx yyy
</pre>
It doesn't matter if the file is the whole protein or just the non-standard residue as the missing parameters should be the same (this is a good check to ensure there are no problems round the corner). The Green Fluorescent Protein (GFP) non-standard residue for example produces the following output:
<pre>
> saveamberparm csy xxx yyy
Checking Unit.
Building topology.
Building atom parameters.
Building bond parameters.
Could not find bond parameter for: CM - HC
Could not find bond parameter for: CC - N*
Could not find bond parameter for: CC - O
Could not find bond parameter for: CC - CM
Could not find bond parameter for: CC - CC
Could not find bond parameter for: CC - N*
Building angle parameters.
Could not find angle parameter: CM - C - OH
Could not find angle parameter: CA - C - CM
Could not find angle parameter: CA - CA - CM
Could not find angle parameter: CA - CM - HC
Could not find angle parameter: CM - CA - CM
Could not find angle parameter: CM - CA - CA
Could not find angle parameter: N3 - CT - H1
Could not find angle parameter: N* - CC - CT
Could not find angle parameter: N* - CT - C
Could not find angle parameter: O - CC - N*
Could not find angle parameter: CC - CC - CM
Could not find angle parameter: CC - N* - CT
Could not find angle parameter: CC - CM - CA
Could not find angle parameter: CC - CM - HC
Could not find angle parameter: CC - CC - N*
Could not find angle parameter: CC - CC - O
Could not find angle parameter: NB - CC - N*
Could not find angle parameter: NB - CC - CM
Could not find angle parameter: NB - CC - CC
Could not find angle parameter: CC - NB - CC
Could not find angle parameter: CC - N* - CC
Could not find angle parameter: CC - N* - CT
Could not find angle parameter: CC - CT - N3
Could not find angle parameter: CC - CT - H1
Building proper torsion parameters.
** No torsion terms for CT-N*-CC-CT
** No torsion terms for N*-CC-CC-CM
** No torsion terms for O-CC-CC-CM
** No torsion terms for O-CC-N*-CT
** No torsion terms for CC-CC-CM-CA
** No torsion terms for CC-CC-CM-HC
** No torsion terms for CC-N*-CC-CT
** No torsion terms for CC-CC-N*-CT
** No torsion terms for NB-CC-N*-CC
** No torsion terms for NB-CC-N*-CT
** No torsion terms for NB-CC-CM-CA
** No torsion terms for NB-CC-CM-HC
** No torsion terms for NB-CC-CC-N*
** No torsion terms for NB-CC-CC-O
** No torsion terms for CC-N*-CC-CC
** No torsion terms for CC-N*-CC-O
Building improper torsion parameters.
total 4 improper torsions applied
Building H-Bond parameters.
Parameter file was not saved.
</pre>
This must then be set up in the Gaussian input file, two lines after the connectivity, in the following style:
=====Bonds=====
{| class="wikitable" border="1"
|-
!'''HrmStr1''': Harmonic stretch I (Amber 1): <math>ForceC(R-R_{eq})^{2}</math>
|-
|HrmStr1 Atom-type1 Atom-type2 ''ForceC'' <math>R_{eq}</math>

''ForceC'' Force constant

<math>R_{eq}</math> Equilibrium bond length
|}

=====Angles=====
{| class="wikitable" border="1"
|-
!'''HrmBnd1''': Harmonic bend (Amber 1): <math>ForceC(\theta-\theta_{eq})^{2}</math>
|-
|HrmBnd1 Atom-type1 Atom-type2 Atom-type3 ''ForceC'' <math>\theta_{eq}</math>

''ForceC'' Force constant (<math>Kcal mol^{-1}rad^{-2}</math>)

<math>\theta_{eq}</math> Equilibrium angle
|}

=====Torsions=====
{| class="wikitable" border="1"
|-
!'''AmbTrs''': Amber torsion (Amber 1): <math>\sum_{i=1}^{4} \frac{Mag_{i}[1+\cos(i\theta - POI(i+4))]}{N_{Paths}}</math>
|-
|AmbTrs Atom-type1 A-type2 A-type3 A-type4 ''PO1 PO2 PO3 PO4'' <math>Mag_{1}</math> <math>Mag_{2}</math> <math>Mag_{3}</math> <math>Mag_{4}</math> <math>N_{Paths}</math>

''PO1–PO4'' Phase offsets for <math>theta</math>: these may be set to 0 or 180: in the former case, they have no effect, in the latter, they have the sole effect of switching the sign of the '+1' coefficient in front of cos.

<math>Mag_{1}</math>-<math>Mag_{4}</math> <math>frac{V}{2}</math> magnitudes

<math>N_{Paths}</math> Number of paths. When zero or less, determined on-the-fly.
|}

Thus for the above example we would obtain a list that looks like:
<pre>
HrmStr1 CM HC
HrmStr1 CC N*
HrmStr1 CC O
HrmStr1 CC CM
HrmStr1 CC CC
HrmBnd1 CM C OH
HrmBnd1 CA C CM
HrmBnd1 HC CM CA
HrmBnd1 CM CA CM
HrmBnd1 CM CA CA
HrmBnd1 H1 CT N3
HrmBnd1 N* CC CT
HrmBnd1 N* CT C
HrmBnd1 O CC N*
HrmBnd1 CC CC CM
HrmBnd1 CC CM CA
HrmBnd1 CC CT N
HrmBnd1 CC CM HC
HrmBnd1 CC CC N*
HrmBnd1 CC CC O
HrmBnd1 NB CC N*
HrmBnd1 NB CC CM
HrmBnd1 NB CC CC
HrmBnd1 CC NB CC
HrmBnd1 CC N* CC
HrmBnd1 CC N* CT
HrmBnd1 CC CT N3
HrmBnd1 CC CT H1
AmbTrs CT N* CC CT
AmbTrs N* CC CC CM
AmbTrs O CC CC CM
AmbTrs O CC N* CT
AmbTrs CC CC CM CA
AmbTrs CC CC CM HC
AmbTrs CC N* CC CT
AmbTrs CC CC N* CT
AmbTrs NB CC N* CC
AmbTrs NB CC N* CT
AmbTrs NB CC CM CA
AmbTrs NB CC CM HC
AmbTrs NB CC CC N*
AmbTrs NB CC CC O
AmbTrs CC N* CC CC
AmbTrs CC N* CC O
</pre>
Now all that remains is to add the values to these parameters. To do this we go back to the files output when we ran the Redtools jobscript. Take the Mol-m1-o1-sm-mol2 file and open it in Gaussview. Change the PDB atom name and AMBER atom types of the .mol2 file in a text editor to those shown in the Gaussview atom list. This should be similar to the file below, obtained for GFP:

[[Image:Mol2_paramget.png|500 px|alt=File for getting AMBER parameters]]

Now we have this file we can obtain the missing parameters from the General AMBER Force Field using:
<pre>
parmchk -i Mol-m1-o1-sm.mol2 -f mol2 -o filename.frcmod
</pre>
where filename can be whatever you chose. The output of this file should now contain all the parameters required for the non-standard residues and the labels for the AMBER atom types should correspond directly to those output by ''saveamberparm''. For stretches and bends the numbers can be simply copied across, however, the torsions are a bit more complicated. An example of a torsion parameter from the .frcmod is shown below:
<pre>
H1-CT-C -O 1 0.800 0.000 -1.000 same as hc-c3-c -o
H1-CT-C -O 1 0.080 180.000 3.000 same as hc-c3-c -o
</pre>
Which must be put in the format

AmbTrs Atom-type1 A-type2 A-type3 A-type4 ''PO1 PO2 PO3 PO4'' <math>Mag_{1}</math> <math>Mag_{2}</math> <math>Mag_{3}</math> <math>Mag_{4}</math> <math>N_{Paths}</math>

Here H1, CT, C and O are the atom types; <math>N_{Paths}</math> is the second column; <math>mag_{i}</math> is the third column; ''POI'' is the fourth column; and the fifth column is the value of i/I. If there is a dash marker, this means that the next row is of the same torsion. The above example would translate then as:
<pre>
AmbTrs H1 CT C O 0 0 180 0 0.8000 0.0000 0.0800 0.0000 1.0
</pre>

We have now determined all the parameters for the AMBER calculation. In order to use them add ''amber=softfirst'' in the route section of the input file. A final point is that in the Gaussian input parameters, the atoms can be specified either way round (e.g. H1 CT C O or O C CT H1). These are equivalent and the input must be checked to ensure that each specification is unique, otherwise an error message will result. This happens even if the values are equal.

A example input file for GFP is shown here: [[Media:1W7S_01_SPE_amber.gjf|Gaussian AMBER input for GFP]]

{| class="wikitable" border="1"
|-
!Known Gaussian and Gaussview labelling problems
|-
|Hydroxyl protons are specified as HO in the AMBER atom types and in the parm96 force field parameters they have zero van der Waals radius. This results in Gaussian showing a warning that charged centres with zero van der Waals radii can collapse into a nearby oppositely charged centre, however, these centres should not have any radii associated with them.
|-
|Carbonyl oxygens are often specified as 'OM' in Gaussview, including any .com files it outputs. This is not recognized in any AMBER parameter sets, however, it is the same as 'O'. Any instances of 'OM' should be changed to 'O'. If this is not done the same problems of charged centres with no van der Waals radii as above occur.
|-
|The backbone nitrogens are often incorrectly labelled as 'N3', however, they should be labelled 'N'.
|-
|Aromatic carbon atoms are labelled 'CH', however, they should be 'CA'.
|}

==Constructing the ONIOM input==
We now have an Gaussian input file which will produce an AMBER calculation of the structure originally specified in the PDB database. This final section details how to progress from this point to an ONIOM input file. In order to check that everything is in order it may be worth running a single point AMBER calculation on the structure. This can be done using ''IOp(4/119=10)'' which will print out the force field parameters so that they can be checked. Providing there are no missing parameters this calculation should complete and an AMBER optimization can be carried out on this structure.

The converged AMBER structure can then be used to make the ONIOM file. A problem here is that there is no partial charge data upon opening the .log file in Gaussview, and the formatted checkpoint file loses all the PDB data. The easiest way to solve this then is to save the .log file as a .pdb file and then to open the .pdb file in Gaussview and save it as a .com file using Calculate→Gaussian Calculation Setup, although this will require re-entering the charges on the non-standard residue. Another option may be to use this [[Media:Onistruct.sh|script]] to add PDB data to the .com file obtained from the .fchk file. Although this is not necessary here, it is useful if you want to create an input with a geometry obtained from a calculation that used geom=check as the .log file loses all PDB data as well.

Whichever option you choose, open the file in Gaussview, specify the high level region using Edit→Edit Layer and then save as a .com file (using Calculate→Gaussian Calculation Setup as otherwise all MM charge data will be lost). The route section you use should look like this for mechanical embedding:
<pre>
#p opt oniom(b3lyp/6-31g(d):amber=softfirst) geom=connectivity
</pre>
or this for electronic embedding:
<pre>
#p opt oniom(b3lyp/6-31g(d):amber=softfirst)=embed geom=connectivity
</pre>

We now have a complete ONIOM input file: [[Media:1W7S_01_OPT_oniom_b3lyp_631gd_amber.gjf|1W7S_01_OPT_oniom_b3lyp_631gd_amber.gjf]]

If you try to run this calculation and get a missing parameter error, this is highly likely to be due to the fact that the parameters involving the link atoms are not present (these parameters could be obtained in the previous steps by using the actual model structure rather that the model without link atoms as done above and the above method should be modified to do this (perhaps using antechamber on a pdb file from Gaussview) in due course). If this happens, check that the missing parameters are not an indication of bad ONIOM partitioning (such as a link-atom replacing an electron-withdrawing group) and add the parameters by hand from the General AMBER Force Field parameter set (available from the AMBER website).

==Summary==
[[Image:Flow_chart_PDB_to_ONIOM.jpg|700 px|Flow-chart summarising procedure for turning PDB structure file to ONIOM input file.]]

Back to [https://wiki.ch.ic.ac.uk/wiki/index.php?title=ONIOM_for_biomolecules ONIOM for biomolecules]

Guide to Creating ONIOM input files for biomolecules

2012-09-05T11:07:37Z

Lmt09: /* Producing a Gaussian Input File */

==Overview==

This guide provides a step by step process to create ONIOM input files for biomolecules from a structure file in the [http://www.pdb.org/pdb/home/home.do Protein Data Bank]. Emphasis will be placed on the use of [https://wiki.ch.ic.ac.uk/wiki/index.php?title=AMBER AMBER] as the low level method and how to obtain parameters for any non-standard residues.

==Creating Standarized .pdb Files==

The first step is to select a .pdb file from the Protein Data Bank that is of high enough resolution to allow atomistic calculations to be produced. The relevant data for determining this is shown on the right hand side under experimental details. The two values to look at are the Resolution[Å] and R-Value, which both should be as low as possible. Having determined a suitable structure, download the suitable test pdb file (usually found in the download files drop-down menu in the top right corner).

In Gaussview select File→Open and choose options. Change the drop-down box "Add Hydrogens:" to '''No''' and, if you wish to remove water molecules, '''check the box''' "Skip Water Molecules." When the file opens up there may be a number of secondary structures present labelled A, B, C etc. In general we require only one so any extras can be removed using Edit→PDB Secondary Structure... and deleting those which are not required. This was then saved as a .pdb file.

{| class="wikitable" border="1"
|-
!The Secondary Structure Editor
|-
|An example of the secondary structure editor is shown below. If we were interested in obtaining structure A only then all that is required is to highlight chains B to D and Edit→Delete→Delete Selected Secondary Structures. The numbered residues such as Helix and Sheet that do not belong to A are automatically removed so if you remove these separately you may end up removing residues from the structure you wish to keep.
[[Image:SSE_PDB_ONIOM_WIKI.png|800 px|alt=PDB Secondary Structure Editor]]
|}

===Residue Names and Protonation States===
Within the .pdb file the fourth column corresponds to the residue name. This name will be used to define the protonation state of the residue, which is currently specified as a default value. In order to determine the protonation state it is possible to use either [http://www.poissonboltzmann.org/pdb2pqr/ PROPKA] or [http://biophysics.cs.vt.edu/ H++]. Once the protonation states have been determined the residue names can be changed to reflect this. (Note: Parameters for non-standard residues calculated later may be included to improve the accuracy).

===Chromophore Structure===
Now a standardized .pdb file of the whole protein has been created the next step is to obtain a .pdb file of the non-standard residue. To do this open the .pdb file we have just saved using a text editor such as vim and remove all lines that are not atoms from the region we intend to include in this residue. It is important here to consider exactly what this consists of here as any problems at this stage are normally not highlighted until much later in the process and will require returning to this point. The region specified here is not the same as that of the ONIOM model region or even the protein chromophore, it is simply so that non-standard residues are defined in the AMBER program. The two important points are that this region must:
# Include the non-standard residue that requires parameterization.
# Is connected to the rest of the protein through standard N or C amino terminations.

The second point may require some elaboration. Some non-standard residues are a modified standard residue, such as that in PYP which is a cystine residue with ''p''-coumaric acid group on the sulphur instead of a thiol. It is tempting to specify the chromophore as just the ''p''-coumaric acid group, however, this causes problems later in defining the parameters for the cystine residue and so the cystine group must also be included in the chromophore region. This joins to the rest of the protein through standard amino acid N and C bonds and so this is all that is needs to be included.

This structure is then saved as a .pdb file and opened in Gaussview. Hydrogens were then added to the residue except where the residue will join to the protein structure. Again be sure of the protonation at this stage as any mistakes will require returning to this point. Check particularly the multiplicity is correct. Save this as a .pdb file and inspect it to ensure that the newly added hydrogens have the same pdb residue name and number as the other atoms, and that their atom numbers follow on and are consistent with connectivity. Also remove any extra TER lines other than the one at the bottom (if there is one). To ensure that this is absolutely correct it may be worth opening this in Gaussview and re-saving it, making sure the correct connectivity is shown.

==Obtaining AMBER Library File of the Chromophore==

We now have two .pdb files, one of the whole protein and one of the non-standard residue region. The next step is to create an AMBER library file of this non-standard residue. Leap, an AMBERTools program, will be used and this requires us to determine three pieces of information for the non-standard residue:
#Connectivity
#AMBER atom types
#Partial charges

Leap can be opened using the command:
<pre>
xleap -s -f /apps/ambertools/amber11/dat/leap/cmd/leaprc.ff03 &
</pre>
This command will use the leap.ff03 set of parameters although any other AMBER parameters could be used depending on the system under study. If this doesn't do anything you probably need to load ambertools:
<pre>
moduleload ambertools
</pre>
Having opened Leap the non-standard residue .pdb file can be loaded using:
<pre>
variable = loadpdb filename
</pre>
where variable is any name you choose and the full pathname must be specified in the filename. Now type:
<pre>
edit variable
</pre>
This brings up a gui where the residue can be visualized. Ensure all atoms are selected and go to Edit→Edit Selected Atoms. This provides a table to be filled with the information specified above. The way to obtain these values will now be explained. A quick sidenote, do not close any Leap x-windows, other than using File→Quit as this will cause the program to crash and any unsaved information to be lost.

====Connectivity====
This is simply achieved by selecting the '''draw''' checkbox in the Leap GUI tool and drawing bonds between the atom centres as desired.

[[Image:Conn_Leap.png|400 px|alt=Drawing bonds with the Leap GUI]]

====AMBER Atom Types====
In order to obtain these open the non-standard residue .pdb file with Gaussview and add methyl groups to the atoms which were previously left with free valences. Save this structure as a .pdb file as we will need it later, however, at this point we only need to go to Edit→Atom List and look at AMBER Type. Copy these across to the Leap table using the PDB Atom Name column to match up Atoms.

====Partial Charges====
This is the most complicated process and requires the use of [[http://q4md-forcefieldtools.org/RED/ R.E.D.-III.4 tools]]. This first uses the modified .pdb file with added methyls to obtain a Gaussian input file using the following command:
<pre>
perl $DIR1/Ante_Red.pl $DIR2/modified_non_standard_residue_file.pdb >> $DIR3/output.txt
</pre>
where $DIR is the relevant pathname. The resulting Gaussian input file can then be run (remember to change memory requirements and checkpoint file locations before submitting). After this has completed the frequency portion was deleted from the log file (this could be removed from the input but is useful for ensuring a minima is obtained) and the log file was copied to Mol_red1.log file in the RED-III directory, ensuring that the filename remains Mol_red1.log. Another file that was output from the above command was a .p2n file. This must be copied to Mol_red1.p2n in the same directory as before, also maintaining Mol_red1.p2n as the filename.

Moving to the RED-III directory now, open Mol_red1.p2n with a text editor and add the following line to exclude the methyl groups from the Partial Charge calculation:
<pre>
REMARK INTRA-MCC 0.0 | 29 30 31 32 33 34 35 36 | R
</pre>
where the numbers correspond to the numbers of the atoms in the methyl groups of the modified non-standard residue. Note that there are two spaces between all the numbers. Below is an example of where it has been placed:

[[Image:p2n_part_char.png|500 px|alt=.p2n file for the calculation of partial charges]]

If necessary change the charge and multiplicity here. Having done this open RED-vIII.4.pl and go to line 4196. Change the variable $DIR to whatever you wish, this is where the output files will be saved to. Create the following jobscript file and run it, although change the directories on line 19 to something useful for you.

<pre>
##################################################################
# REDTOOLS JOBSCRIPT #
# CREATED 08/07/10 #
# LAST MODIFIED 08/07/10 #
# LEE THOMPSON #
##################################################################

#PBS -l ncpus=1
#PBS -l mem=1000mb
#PBS -l walltime=04:00:00
#PBS -joe

module load ambertools
module load gaussian
export GAUSS_SCRDIR=$TMPDIR
echo $GAUSS_SCRDIR
cd $(echo $PBS_O_WORKDIR)
pwd
perl /home/lmt09/SOFTWARE/RED-III.4-Tools-Files/RED-vIII.4.pl > /home/lmt09/PHD_Y2/PYP/1NWZ/PROTONATED/ONIOM/RED_out.log

</pre>

Going into the new directory open Mol-m1-o1-sm.mol2 which contains the partial charges that we seek in the final column. To copy these to the Leap table requires a bit of detective work to match up the atoms. This can be done by opening up the Gaussian log file Mol_red1.log in Gaussview which is labelled in the same order as the .mol2 file with the partial charges. The Gaussview atoms and the atoms in the Leap GUI can then be matched by their positions. These atoms can then be matched to the Leap table by displaying atom names on the Leap GUI using Display→Names. This is also a good time to check consistency of atom types again as if they are different it will cause problems identifying parameters later on. It is also worth checking that the charges sum to an integer value and that '''you have typed them in correctly'''.

----

Having filled in the Leap table go to File→Save and Quit, and then exit the GUI using File→Close. Back at the command line prompt, save the library file using:

<pre>
saveoff variable filename
</pre>
where the variable is the same as used before and filename includes the full pathname. Now exit Leap and go to the .lib file that we have just created. In order for this to be recognised the filename must be uppercase and three or four letters long (although I have not tried to see otherwise). In order to achieve this move it from variable.lib to VAR.lib, where, VAR is a capitalized three letter word of your choice. Now open the file in vi and type:
<pre>
:%s/variable/VAR/g
</pre>
which changes all instances of variable to VAR. We have now created our library file for the non-standard residue.

{| class="wikitable" border="1"
|-
!Justification of Partial Charge Model
|-
|The determination of partial charges is important for the successful use of force field methods, yet the concept of a partial charge is somewhat ambiguous, with several different methods for their determination (see Cramer, C.J., ''Essentials of Computational Chemistry, p309'' for an introduction. The partial charges we use are computed using the restrained ESP method (Cornell ''et al'', ''Journal of the American Chemical Society'', '''1995''', 117, 19, 5179-5197). This is an extension of the ESP method which determines partial charge ''q'' on atom ''k'' by minimizing the difference between:
<math>V_{ESP}(\mathbf{r}) = \sum_{k}^{nuclei}\frac{q_{k}}{\vert \mathbf{r}-\mathbf{r}_{k}\vert}</math><br\>
and the Molecular Electrostatic Potential (MEP):<br\>
<math>V_{MEP}(\mathbf{r}) = \sum_{k}^{nuclei}\frac{Z_{k}}{\vert \mathbf{r}-\mathbf{r}_{k}\vert}-\int \! \Psi (\mathbf{r'}) \frac{1}{\vert \mathbf{r}-\mathbf{r'}\vert} \Psi (\mathbf{r'}) \, \mathrm{d}\mathbf{r'}</math><br\>
for all positions '''r'''. This is computed from a number of points spaced evenly around the Connolly surface of the molecule. ESP is dependent on conformation, however, causing hydrogens in a methyl group for example to have different partial charges. As these are all freely rotating in practice the same partial charges may used for each hydrogen and this is the extension that RESP applies to the ESP method (Bayly ''et al'', ''The Journal of Physical Chemistry'', '''1993''', 97, 40, 10269-10280).
|-
|The main purpose for using this is that AMBER uses RESP for its parm96 (Cornell) parameter set which is the same as that used by Gaussian (derived from HF/6-31G*). Reasons for its use in this force field are that it has been shown to be useful for modeling inter-molecular interactions at short to long range, is convergent with respect to the size of basis set used, resolves to an extent the problems of atoms which do not contribute the Connolly surface and so are ill-defined by the method, as well as having the original advantages of ESP over methods such as Mulliken and Löwdin charges.
|-
|RED (RESP and ESP charge Derive) tools is a series of perl scripts which generate a Gaussian input file which can be run and from which the partial charges derived (Dupradeau ''et al'', ''Physical chemistry chemical physics: PCCP'', '''2010''', 12, 28, 7821-39).
|}

==Producing a Gaussian Input File==

'''NOTE: It is now recomemnded to add hydrogens using PROPKA or H++ prior to loading into LEAP'''

Having constructed the library file of the non-standard residue we must now construct a .com or .gjf file to run in Gaussian. Initially this will simply be an AMBER calculation, the output of which will be used to determine if we have all the correct parameters and as a starting geometry for the ONIOM calculations. The first step is to reopen Leap using the same command as before. now load in the AMBER library file for the non standard residue using the command:
<pre>
loadoff $DIR/VAR.lib
</pre>
where again $DIR represents the pathname of the file. Typing ''list'' in Leap will display all the library files that have been loaded of which VAR should be one of them. The next stage is to load the .pdb file of the protein that we obtained from Gaussview previously using the command for loading .pdb files shown previously. This should add hydrogens to the structure in accordance to the library files and perhaps a terminal oxygen although never any other heavy atom (this is displayed at the command line). Opening the Leap GUI of the whole protein should reveal the non-standard residue highlighted in the full protein environment. The connection between the protein and the residue must now be determined. There are a number of ways to do this including at the command line but the most successful method so far is to simply draw the bonds in the Leap GUI.

[[Image:Pro_struct.png|500 px|alt=Protein structure in Leap GUI]]

Having drawn the connectivity, go to Unit→Calculate Net Charge to obtain the charge of the protein, which should be an integer. Close the Leap GUI and save the .pdb using the following command:
<pre>
savepdb variable $DIR/filename.pdb
</pre>
We will need the information we have loaded in Leap later so do not close it but for now, look at the .pdb file in the text editor and ensure that there is a terminal oxygen labelled <OXT> at the bottom of the file. If there is not insert ''ATOM 3608 OXT HIE 228 -1.012 21.725 100.791 1.00 0.00'' in the correct place, although the cartesian coordinates, PDB residue name and number and atom number will be different from this example. Also bear in mind that if there are any waters below this then their atom numbers will need changing (use ''grep "WAT" filename1.pdb | awk '{ X=$2; Y=X+1; print "s/"X," "$3" /"Y," "$3" /" }' > sedscript | sed -f sedscript < filename1.pdb > filename2.pdb'' for this). Also check in the .pdb file that the atoms around the chromophore all are of the same residue and do not have differing residue numbers. If this is not the case then it means that the non-standard residue .pdb file was mixed up and you must return to that stage.

Now open the .pdb file in Gaussview and go to Edit→Atom List. Scan through this to ensure that all MM partial charges are present for all atoms other than those in the chromophore residue. If there are any that are not it is because the residue connectivity is wrong so use the bond specification tool to correct this in Gaussview and you should see the MM charges appear as soon as you correct the problem (Hint: The problem atoms will have undefined AMBER atom types (shown as ?) so look at connectivity around these atoms). The MM partial charges can be copied directly from Mol-m1-o1-sm.mol2 into the Gaussview atom list now, although I prefer a second option which I shall explain when I come to it. Now save this as a .com/.gjf file although, because of a bug in Gaussview which causes patial charges to be missing from the input file, you '''must''' save this using Calculate→Gaussian Calculation Setup, chose an AMBER calculation and insert the charge determined earlier and the multiplicity and submit. Select yes when prompted to save the file and then cancel the file execution. You should now have a Gaussian input file in your directory. If you have not inserted the MM partial charges previously copy them from Mol-m1-o1-sm.mol2 and paste them between the AMBER atom type and the PDB information in the Gaussian input file. All that we now require for a complete AMBER calculation is the AMBER parameters for the non-standard residue.

==Getting Non-Standard AMBER Parameters==
If we were to run the Gaussian input file as produced above we would get an error message indicating missing AMBER parameters. Gaussian uses parm96 by default and if any stretches, bends or torsions are present in the non-standard residue but not in the forcefield, then an error message is obtained. In order to obtain them we first need to know what parameters are missing. This can be achieved using:
<pre>
saveamberparm variable xxx yyy
</pre>
It doesn't matter if the file is the whole protein or just the non-standard residue as the missing parameters should be the same (this is a good check to ensure there are no problems round the corner). The Green Fluorescent Protein (GFP) non-standard residue for example produces the following output:
<pre>
> saveamberparm csy xxx yyy
Checking Unit.
Building topology.
Building atom parameters.
Building bond parameters.
Could not find bond parameter for: CM - HC
Could not find bond parameter for: CC - N*
Could not find bond parameter for: CC - O
Could not find bond parameter for: CC - CM
Could not find bond parameter for: CC - CC
Could not find bond parameter for: CC - N*
Building angle parameters.
Could not find angle parameter: CM - C - OH
Could not find angle parameter: CA - C - CM
Could not find angle parameter: CA - CA - CM
Could not find angle parameter: CA - CM - HC
Could not find angle parameter: CM - CA - CM
Could not find angle parameter: CM - CA - CA
Could not find angle parameter: N3 - CT - H1
Could not find angle parameter: N* - CC - CT
Could not find angle parameter: N* - CT - C
Could not find angle parameter: O - CC - N*
Could not find angle parameter: CC - CC - CM
Could not find angle parameter: CC - N* - CT
Could not find angle parameter: CC - CM - CA
Could not find angle parameter: CC - CM - HC
Could not find angle parameter: CC - CC - N*
Could not find angle parameter: CC - CC - O
Could not find angle parameter: NB - CC - N*
Could not find angle parameter: NB - CC - CM
Could not find angle parameter: NB - CC - CC
Could not find angle parameter: CC - NB - CC
Could not find angle parameter: CC - N* - CC
Could not find angle parameter: CC - N* - CT
Could not find angle parameter: CC - CT - N3
Could not find angle parameter: CC - CT - H1
Building proper torsion parameters.
** No torsion terms for CT-N*-CC-CT
** No torsion terms for N*-CC-CC-CM
** No torsion terms for O-CC-CC-CM
** No torsion terms for O-CC-N*-CT
** No torsion terms for CC-CC-CM-CA
** No torsion terms for CC-CC-CM-HC
** No torsion terms for CC-N*-CC-CT
** No torsion terms for CC-CC-N*-CT
** No torsion terms for NB-CC-N*-CC
** No torsion terms for NB-CC-N*-CT
** No torsion terms for NB-CC-CM-CA
** No torsion terms for NB-CC-CM-HC
** No torsion terms for NB-CC-CC-N*
** No torsion terms for NB-CC-CC-O
** No torsion terms for CC-N*-CC-CC
** No torsion terms for CC-N*-CC-O
Building improper torsion parameters.
total 4 improper torsions applied
Building H-Bond parameters.
Parameter file was not saved.
</pre>
This must then be set up in the Gaussian input file, two lines after the connectivity, in the following style:
=====Bonds=====
{| class="wikitable" border="1"
|-
!'''HrmStr1''': Harmonic stretch I (Amber 1): <math>ForceC(R-R_{eq})^{2}</math>
|-
|HrmStr1 Atom-type1 Atom-type2 ''ForceC'' <math>R_{eq}</math>

''ForceC'' Force constant

<math>R_{eq}</math> Equilibrium bond length
|}

=====Angles=====
{| class="wikitable" border="1"
|-
!'''HrmBnd1''': Harmonic bend (Amber 1): <math>ForceC(\theta-\theta_{eq})^{2}</math>
|-
|HrmBnd1 Atom-type1 Atom-type2 Atom-type3 ''ForceC'' <math>\theta_{eq}</math>

''ForceC'' Force constant (<math>Kcal mol^{-1}rad^{-2}</math>)

<math>\theta_{eq}</math> Equilibrium angle
|}

=====Torsions=====
{| class="wikitable" border="1"
|-
!'''AmbTrs''': Amber torsion (Amber 1): <math>\sum_{i=1}^{4} \frac{Mag_{i}[1+\cos(i\theta - POI(i+4))]}{N_{Paths}}</math>
|-
|AmbTrs Atom-type1 A-type2 A-type3 A-type4 ''PO1 PO2 PO3 PO4'' <math>Mag_{1}</math> <math>Mag_{2}</math> <math>Mag_{3}</math> <math>Mag_{4}</math> <math>N_{Paths}</math>

''PO1–PO4'' Phase offsets for <math>theta</math>: these may be set to 0 or 180: in the former case, they have no effect, in the latter, they have the sole effect of switching the sign of the '+1' coefficient in front of cos.

<math>Mag_{1}</math>-<math>Mag_{4}</math> <math>frac{V}{2}</math> magnitudes

<math>N_{Paths}</math> Number of paths. When zero or less, determined on-the-fly.
|}

Thus for the above example we would obtain a list that looks like:
<pre>
HrmStr1 CM HC
HrmStr1 CC N*
HrmStr1 CC O
HrmStr1 CC CM
HrmStr1 CC CC
HrmBnd1 CM C OH
HrmBnd1 CA C CM
HrmBnd1 HC CM CA
HrmBnd1 CM CA CM
HrmBnd1 CM CA CA
HrmBnd1 H1 CT N3
HrmBnd1 N* CC CT
HrmBnd1 N* CT C
HrmBnd1 O CC N*
HrmBnd1 CC CC CM
HrmBnd1 CC CM CA
HrmBnd1 CC CT N
HrmBnd1 CC CM HC
HrmBnd1 CC CC N*
HrmBnd1 CC CC O
HrmBnd1 NB CC N*
HrmBnd1 NB CC CM
HrmBnd1 NB CC CC
HrmBnd1 CC NB CC
HrmBnd1 CC N* CC
HrmBnd1 CC N* CT
HrmBnd1 CC CT N3
HrmBnd1 CC CT H1
AmbTrs CT N* CC CT
AmbTrs N* CC CC CM
AmbTrs O CC CC CM
AmbTrs O CC N* CT
AmbTrs CC CC CM CA
AmbTrs CC CC CM HC
AmbTrs CC N* CC CT
AmbTrs CC CC N* CT
AmbTrs NB CC N* CC
AmbTrs NB CC N* CT
AmbTrs NB CC CM CA
AmbTrs NB CC CM HC
AmbTrs NB CC CC N*
AmbTrs NB CC CC O
AmbTrs CC N* CC CC
AmbTrs CC N* CC O
</pre>
Now all that remains is to add the values to these parameters. To do this we go back to the files output when we ran the Redtools jobscript. Take the Mol-m1-o1-sm-mol2 file and open it in Gaussview. Change the PDB atom name and AMBER atom types of the .mol2 file in a text editor to those shown in the Gaussview atom list. This should be similar to the file below, obtained for GFP:

[[Image:Mol2_paramget.png|500 px|alt=File for getting AMBER parameters]]

Now we have this file we can obtain the missing parameters from the General AMBER Force Field using:
<pre>
parmchk -i Mol-m1-o1-sm.mol2 -f mol2 -o filename.frcmod
</pre>
where filename can be whatever you chose. The output of this file should now contain all the parameters required for the non-standard residues and the labels for the AMBER atom types should correspond directly to those output by ''saveamberparm''. For stretches and bends the numbers can be simply copied across, however, the torsions are a bit more complicated. An example of a torsion parameter from the .frcmod is shown below:
<pre>
H1-CT-C -O 1 0.800 0.000 -1.000 same as hc-c3-c -o
H1-CT-C -O 1 0.080 180.000 3.000 same as hc-c3-c -o
</pre>
Which must be put in the format

AmbTrs Atom-type1 A-type2 A-type3 A-type4 ''PO1 PO2 PO3 PO4'' <math>Mag_{1}</math> <math>Mag_{2}</math> <math>Mag_{3}</math> <math>Mag_{4}</math> <math>N_{Paths}</math>

Here H1, CT, C and O are the atom types; <math>N_{Paths}</math> is the second column; <math>mag_{i}</math> is the third column; ''POI'' is the fourth column; and the fifth column is the value of i/I. If there is a dash marker, this means that the next row is of the same torsion. The above example would translate then as:
<pre>
AmbTrs H1 CT C O 0 0 180 0 0.8000 0.0000 0.0800 0.0000 1.0
</pre>

We have now determined all the parameters for the AMBER calculation. In order to use them add ''amber=softfirst'' in the route section of the input file. A final point is that in the Gaussian input parameters, the atoms can be specified either way round (e.g. H1 CT C O or O C CT H1). These are equivalent and the input must be checked to ensure that each specification is unique, otherwise an error message will result. This happens even if the values are equal.

A example input file for GFP is shown here: [[Media:1W7S_01_SPE_amber.gjf|Gaussian AMBER input for GFP]]

{| class="wikitable" border="1"
|-
!Known Gaussian and Gaussview labelling problems
|-
|Hydroxyl protons are specified as HO in the AMBER atom types and in the parm96 force field parameters they have zero van der Waals radius. This results in Gaussian showing a warning that charged centres with zero van der Waals radii can collapse into a nearby oppositely charged centre, however, these centres should not have any radii associated with them.
|-
|Carbonyl oxygens are often specified as 'OM' in Gaussview, including any .com files it outputs. This is not recognized in any AMBER parameter sets, however, it is the same as 'O'. Any instances of 'OM' should be changed to 'O'. If this is not done the same problems of charged centres with no van der Waals radii as above occur.
|-
|The backbone nitrogens are often incorrectly labelled as 'N3', however, they should be labelled 'N'.
|-
|Aromatic carbon atoms are labelled 'CH', however, they should be 'CA'.
|}

==Constructing the ONIOM input==
We now have an Gaussian input file which will produce an AMBER calculation of the structure originally specified in the PDB database. This final section details how to progress from this point to an ONIOM input file. In order to check that everything is in order it may be worth running a single point AMBER calculation on the structure. This can be done using ''IOp(4/119=10)'' which will print out the force field parameters so that they can be checked. Providing there are no missing parameters this calculation should complete and an AMBER optimization can be carried out on this structure.

The converged AMBER structure can then be used to make the ONIOM file. A problem here is that there is no partial charge data upon opening the .log file in Gaussview, and the formatted checkpoint file loses all the PDB data. The easiest way to solve this then is to save the .log file as a .pdb file and then to open the .pdb file in Gaussview and save it as a .com file using Calculate→Gaussian Calculation Setup, although this will require re-entering the charges on the non-standard residue. Another option may be to use this [[Media:Onistruct.sh|script]] to add PDB data to the .com file obtained from the .fchk file. Although this is not necessary here, it is useful if you want to create an input with a geometry obtained from a calculation that used geom=check as the .log file loses all PDB data as well.

Whichever option you choose, open the file in Gaussview, specify the high level region using Edit→Edit Layer and then save as a .com file (using Calculate→Gaussian Calculation Setup as otherwise all MM charge data will be lost). The route section you use should look like this for mechanical embedding:
<pre>
#p opt oniom(b3lyp/6-31g(d):amber=softfirst) geom=connectivity
</pre>
or this for electronic embedding:
<pre>
#p opt oniom(b3lyp/6-31g(d):amber=softfirst)=embed geom=connectivity
</pre>

We now have a complete ONIOM input file: [[Media:1W7S_01_OPT_oniom_b3lyp_631gd_amber.gjf|1W7S_01_OPT_oniom_b3lyp_631gd_amber.gjf]]

If you try to run this calculation and get a missing parameter error, this is highly likely to be due to the fact that the parameters involving the link atoms are not present (these parameters could be obtained in the previous steps by using the actual model structure rather that the model without link atoms as done above and the above method should be modified to do this (perhaps using antechamber on a pdb file from Gaussview) in due course). If this happens, check that the missing parameters are not an indication of bad ONIOM partitioning (such as a link-atom replacing an electron-withdrawing group) and add the parameters by hand from the General AMBER Force Field parameter set (available from the AMBER website).

==Summary==
[[Image:Flow_chart_PDB_to_ONIOM.jpg|700 px|Flow-chart summarising procedure for turning PDB structure file to ONIOM input file.]]

Back to [https://wiki.ch.ic.ac.uk/wiki/index.php?title=ONIOM_for_biomolecules ONIOM for biomolecules]

Guide to Creating ONIOM input files for biomolecules

2012-09-05T11:06:03Z

Lmt09: /* Producing a Gaussian Input File */

==Overview==

This guide provides a step by step process to create ONIOM input files for biomolecules from a structure file in the [http://www.pdb.org/pdb/home/home.do Protein Data Bank]. Emphasis will be placed on the use of [https://wiki.ch.ic.ac.uk/wiki/index.php?title=AMBER AMBER] as the low level method and how to obtain parameters for any non-standard residues.

==Creating Standarized .pdb Files==

The first step is to select a .pdb file from the Protein Data Bank that is of high enough resolution to allow atomistic calculations to be produced. The relevant data for determining this is shown on the right hand side under experimental details. The two values to look at are the Resolution[Å] and R-Value, which both should be as low as possible. Having determined a suitable structure, download the suitable test pdb file (usually found in the download files drop-down menu in the top right corner).

In Gaussview select File→Open and choose options. Change the drop-down box "Add Hydrogens:" to '''No''' and, if you wish to remove water molecules, '''check the box''' "Skip Water Molecules." When the file opens up there may be a number of secondary structures present labelled A, B, C etc. In general we require only one so any extras can be removed using Edit→PDB Secondary Structure... and deleting those which are not required. This was then saved as a .pdb file.

{| class="wikitable" border="1"
|-
!The Secondary Structure Editor
|-
|An example of the secondary structure editor is shown below. If we were interested in obtaining structure A only then all that is required is to highlight chains B to D and Edit→Delete→Delete Selected Secondary Structures. The numbered residues such as Helix and Sheet that do not belong to A are automatically removed so if you remove these separately you may end up removing residues from the structure you wish to keep.
[[Image:SSE_PDB_ONIOM_WIKI.png|800 px|alt=PDB Secondary Structure Editor]]
|}

===Residue Names and Protonation States===
Within the .pdb file the fourth column corresponds to the residue name. This name will be used to define the protonation state of the residue, which is currently specified as a default value. In order to determine the protonation state it is possible to use either [http://www.poissonboltzmann.org/pdb2pqr/ PROPKA] or [http://biophysics.cs.vt.edu/ H++]. Once the protonation states have been determined the residue names can be changed to reflect this. (Note: Parameters for non-standard residues calculated later may be included to improve the accuracy).

===Chromophore Structure===
Now a standardized .pdb file of the whole protein has been created the next step is to obtain a .pdb file of the non-standard residue. To do this open the .pdb file we have just saved using a text editor such as vim and remove all lines that are not atoms from the region we intend to include in this residue. It is important here to consider exactly what this consists of here as any problems at this stage are normally not highlighted until much later in the process and will require returning to this point. The region specified here is not the same as that of the ONIOM model region or even the protein chromophore, it is simply so that non-standard residues are defined in the AMBER program. The two important points are that this region must:
# Include the non-standard residue that requires parameterization.
# Is connected to the rest of the protein through standard N or C amino terminations.

The second point may require some elaboration. Some non-standard residues are a modified standard residue, such as that in PYP which is a cystine residue with ''p''-coumaric acid group on the sulphur instead of a thiol. It is tempting to specify the chromophore as just the ''p''-coumaric acid group, however, this causes problems later in defining the parameters for the cystine residue and so the cystine group must also be included in the chromophore region. This joins to the rest of the protein through standard amino acid N and C bonds and so this is all that is needs to be included.

This structure is then saved as a .pdb file and opened in Gaussview. Hydrogens were then added to the residue except where the residue will join to the protein structure. Again be sure of the protonation at this stage as any mistakes will require returning to this point. Check particularly the multiplicity is correct. Save this as a .pdb file and inspect it to ensure that the newly added hydrogens have the same pdb residue name and number as the other atoms, and that their atom numbers follow on and are consistent with connectivity. Also remove any extra TER lines other than the one at the bottom (if there is one). To ensure that this is absolutely correct it may be worth opening this in Gaussview and re-saving it, making sure the correct connectivity is shown.

==Obtaining AMBER Library File of the Chromophore==

We now have two .pdb files, one of the whole protein and one of the non-standard residue region. The next step is to create an AMBER library file of this non-standard residue. Leap, an AMBERTools program, will be used and this requires us to determine three pieces of information for the non-standard residue:
#Connectivity
#AMBER atom types
#Partial charges

Leap can be opened using the command:
<pre>
xleap -s -f /apps/ambertools/amber11/dat/leap/cmd/leaprc.ff03 &
</pre>
This command will use the leap.ff03 set of parameters although any other AMBER parameters could be used depending on the system under study. If this doesn't do anything you probably need to load ambertools:
<pre>
moduleload ambertools
</pre>
Having opened Leap the non-standard residue .pdb file can be loaded using:
<pre>
variable = loadpdb filename
</pre>
where variable is any name you choose and the full pathname must be specified in the filename. Now type:
<pre>
edit variable
</pre>
This brings up a gui where the residue can be visualized. Ensure all atoms are selected and go to Edit→Edit Selected Atoms. This provides a table to be filled with the information specified above. The way to obtain these values will now be explained. A quick sidenote, do not close any Leap x-windows, other than using File→Quit as this will cause the program to crash and any unsaved information to be lost.

====Connectivity====
This is simply achieved by selecting the '''draw''' checkbox in the Leap GUI tool and drawing bonds between the atom centres as desired.

[[Image:Conn_Leap.png|400 px|alt=Drawing bonds with the Leap GUI]]

====AMBER Atom Types====
In order to obtain these open the non-standard residue .pdb file with Gaussview and add methyl groups to the atoms which were previously left with free valences. Save this structure as a .pdb file as we will need it later, however, at this point we only need to go to Edit→Atom List and look at AMBER Type. Copy these across to the Leap table using the PDB Atom Name column to match up Atoms.

====Partial Charges====
This is the most complicated process and requires the use of [[http://q4md-forcefieldtools.org/RED/ R.E.D.-III.4 tools]]. This first uses the modified .pdb file with added methyls to obtain a Gaussian input file using the following command:
<pre>
perl $DIR1/Ante_Red.pl $DIR2/modified_non_standard_residue_file.pdb >> $DIR3/output.txt
</pre>
where $DIR is the relevant pathname. The resulting Gaussian input file can then be run (remember to change memory requirements and checkpoint file locations before submitting). After this has completed the frequency portion was deleted from the log file (this could be removed from the input but is useful for ensuring a minima is obtained) and the log file was copied to Mol_red1.log file in the RED-III directory, ensuring that the filename remains Mol_red1.log. Another file that was output from the above command was a .p2n file. This must be copied to Mol_red1.p2n in the same directory as before, also maintaining Mol_red1.p2n as the filename.

Moving to the RED-III directory now, open Mol_red1.p2n with a text editor and add the following line to exclude the methyl groups from the Partial Charge calculation:
<pre>
REMARK INTRA-MCC 0.0 | 29 30 31 32 33 34 35 36 | R
</pre>
where the numbers correspond to the numbers of the atoms in the methyl groups of the modified non-standard residue. Note that there are two spaces between all the numbers. Below is an example of where it has been placed:

[[Image:p2n_part_char.png|500 px|alt=.p2n file for the calculation of partial charges]]

If necessary change the charge and multiplicity here. Having done this open RED-vIII.4.pl and go to line 4196. Change the variable $DIR to whatever you wish, this is where the output files will be saved to. Create the following jobscript file and run it, although change the directories on line 19 to something useful for you.

<pre>
##################################################################
# REDTOOLS JOBSCRIPT #
# CREATED 08/07/10 #
# LAST MODIFIED 08/07/10 #
# LEE THOMPSON #
##################################################################

#PBS -l ncpus=1
#PBS -l mem=1000mb
#PBS -l walltime=04:00:00
#PBS -joe

module load ambertools
module load gaussian
export GAUSS_SCRDIR=$TMPDIR
echo $GAUSS_SCRDIR
cd $(echo $PBS_O_WORKDIR)
pwd
perl /home/lmt09/SOFTWARE/RED-III.4-Tools-Files/RED-vIII.4.pl > /home/lmt09/PHD_Y2/PYP/1NWZ/PROTONATED/ONIOM/RED_out.log

</pre>

Going into the new directory open Mol-m1-o1-sm.mol2 which contains the partial charges that we seek in the final column. To copy these to the Leap table requires a bit of detective work to match up the atoms. This can be done by opening up the Gaussian log file Mol_red1.log in Gaussview which is labelled in the same order as the .mol2 file with the partial charges. The Gaussview atoms and the atoms in the Leap GUI can then be matched by their positions. These atoms can then be matched to the Leap table by displaying atom names on the Leap GUI using Display→Names. This is also a good time to check consistency of atom types again as if they are different it will cause problems identifying parameters later on. It is also worth checking that the charges sum to an integer value and that '''you have typed them in correctly'''.

----

Having filled in the Leap table go to File→Save and Quit, and then exit the GUI using File→Close. Back at the command line prompt, save the library file using:

<pre>
saveoff variable filename
</pre>
where the variable is the same as used before and filename includes the full pathname. Now exit Leap and go to the .lib file that we have just created. In order for this to be recognised the filename must be uppercase and three or four letters long (although I have not tried to see otherwise). In order to achieve this move it from variable.lib to VAR.lib, where, VAR is a capitalized three letter word of your choice. Now open the file in vi and type:
<pre>
:%s/variable/VAR/g
</pre>
which changes all instances of variable to VAR. We have now created our library file for the non-standard residue.

{| class="wikitable" border="1"
|-
!Justification of Partial Charge Model
|-
|The determination of partial charges is important for the successful use of force field methods, yet the concept of a partial charge is somewhat ambiguous, with several different methods for their determination (see Cramer, C.J., ''Essentials of Computational Chemistry, p309'' for an introduction. The partial charges we use are computed using the restrained ESP method (Cornell ''et al'', ''Journal of the American Chemical Society'', '''1995''', 117, 19, 5179-5197). This is an extension of the ESP method which determines partial charge ''q'' on atom ''k'' by minimizing the difference between:
<math>V_{ESP}(\mathbf{r}) = \sum_{k}^{nuclei}\frac{q_{k}}{\vert \mathbf{r}-\mathbf{r}_{k}\vert}</math><br\>
and the Molecular Electrostatic Potential (MEP):<br\>
<math>V_{MEP}(\mathbf{r}) = \sum_{k}^{nuclei}\frac{Z_{k}}{\vert \mathbf{r}-\mathbf{r}_{k}\vert}-\int \! \Psi (\mathbf{r'}) \frac{1}{\vert \mathbf{r}-\mathbf{r'}\vert} \Psi (\mathbf{r'}) \, \mathrm{d}\mathbf{r'}</math><br\>
for all positions '''r'''. This is computed from a number of points spaced evenly around the Connolly surface of the molecule. ESP is dependent on conformation, however, causing hydrogens in a methyl group for example to have different partial charges. As these are all freely rotating in practice the same partial charges may used for each hydrogen and this is the extension that RESP applies to the ESP method (Bayly ''et al'', ''The Journal of Physical Chemistry'', '''1993''', 97, 40, 10269-10280).
|-
|The main purpose for using this is that AMBER uses RESP for its parm96 (Cornell) parameter set which is the same as that used by Gaussian (derived from HF/6-31G*). Reasons for its use in this force field are that it has been shown to be useful for modeling inter-molecular interactions at short to long range, is convergent with respect to the size of basis set used, resolves to an extent the problems of atoms which do not contribute the Connolly surface and so are ill-defined by the method, as well as having the original advantages of ESP over methods such as Mulliken and Löwdin charges.
|-
|RED (RESP and ESP charge Derive) tools is a series of perl scripts which generate a Gaussian input file which can be run and from which the partial charges derived (Dupradeau ''et al'', ''Physical chemistry chemical physics: PCCP'', '''2010''', 12, 28, 7821-39).
|}

==Producing a Gaussian Input File==

''NOTE: It is now reccomended to add hydrogens using PROPKA or H++ prior to loading into LEAP''

Having constructed the library file of the non-standard residue we must now construct a .com or .gjf file to run in Gaussian. Initially this will simply be an AMBER calculation, the output of which will be used to determine if we have all the correct parameters and as a starting geometry for the ONIOM calculations. The first step is to reopen Leap using the same command as before. now load in the AMBER library file for the non standard residue using the command:
<pre>
loadoff $DIR/VAR.lib
</pre>
where again $DIR represents the pathname of the file. Typing ''list'' in Leap will display all the library files that have been loaded of which VAR should be one of them. The next stage is to load the .pdb file of the protein that we obtained from Gaussview previously using the command for loading .pdb files shown previously. This should add hydrogens to the structure in accordance to the library files and perhaps a terminal oxygen although never any other heavy atom (this is displayed at the command line). Opening the Leap GUI of the whole protein should reveal the non-standard residue highlighted in the full protein environment. The connection between the protein and the residue must now be determined. There are a number of ways to do this including at the command line but the most successful method so far is to simply draw the bonds in the Leap GUI.

[[Image:Pro_struct.png|500 px|alt=Protein structure in Leap GUI]]

Having drawn the connectivity, go to Unit→Calculate Net Charge to obtain the charge of the protein, which should be an integer. Close the Leap GUI and save the .pdb using the following command:
<pre>
savepdb variable $DIR/filename.pdb
</pre>
We will need the information we have loaded in Leap later so do not close it but for now, look at the .pdb file in the text editor and ensure that there is a terminal oxygen labelled <OXT> at the bottom of the file. If there is not insert ''ATOM 3608 OXT HIE 228 -1.012 21.725 100.791 1.00 0.00'' in the correct place, although the cartesian coordinates, PDB residue name and number and atom number will be different from this example. Also bear in mind that if there are any waters below this then their atom numbers will need changing (use ''grep "WAT" filename1.pdb | awk '{ X=$2; Y=X+1; print "s/"X," "$3" /"Y," "$3" /" }' > sedscript | sed -f sedscript < filename1.pdb > filename2.pdb'' for this). Also check in the .pdb file that the atoms around the chromophore all are of the same residue and do not have differing residue numbers. If this is not the case then it means that the non-standard residue .pdb file was mixed up and you must return to that stage.

Now open the .pdb file in Gaussview and go to Edit→Atom List. Scan through this to ensure that all MM partial charges are present for all atoms other than those in the chromophore residue. If there are any that are not it is because the residue connectivity is wrong so use the bond specification tool to correct this in Gaussview and you should see the MM charges appear as soon as you correct the problem (Hint: The problem atoms will have undefined AMBER atom types (shown as ?) so look at connectivity around these atoms). The MM partial charges can be copied directly from Mol-m1-o1-sm.mol2 into the Gaussview atom list now, although I prefer a second option which I shall explain when I come to it. Now save this as a .com/.gjf file although, because of a bug in Gaussview which causes patial charges to be missing from the input file, you '''must''' save this using Calculate→Gaussian Calculation Setup, chose an AMBER calculation and insert the charge determined earlier and the multiplicity and submit. Select yes when prompted to save the file and then cancel the file execution. You should now have a Gaussian input file in your directory. If you have not inserted the MM partial charges previously copy them from Mol-m1-o1-sm.mol2 and paste them between the AMBER atom type and the PDB information in the Gaussian input file. All that we now require for a complete AMBER calculation is the AMBER parameters for the non-standard residue.

==Getting Non-Standard AMBER Parameters==
If we were to run the Gaussian input file as produced above we would get an error message indicating missing AMBER parameters. Gaussian uses parm96 by default and if any stretches, bends or torsions are present in the non-standard residue but not in the forcefield, then an error message is obtained. In order to obtain them we first need to know what parameters are missing. This can be achieved using:
<pre>
saveamberparm variable xxx yyy
</pre>
It doesn't matter if the file is the whole protein or just the non-standard residue as the missing parameters should be the same (this is a good check to ensure there are no problems round the corner). The Green Fluorescent Protein (GFP) non-standard residue for example produces the following output:
<pre>
> saveamberparm csy xxx yyy
Checking Unit.
Building topology.
Building atom parameters.
Building bond parameters.
Could not find bond parameter for: CM - HC
Could not find bond parameter for: CC - N*
Could not find bond parameter for: CC - O
Could not find bond parameter for: CC - CM
Could not find bond parameter for: CC - CC
Could not find bond parameter for: CC - N*
Building angle parameters.
Could not find angle parameter: CM - C - OH
Could not find angle parameter: CA - C - CM
Could not find angle parameter: CA - CA - CM
Could not find angle parameter: CA - CM - HC
Could not find angle parameter: CM - CA - CM
Could not find angle parameter: CM - CA - CA
Could not find angle parameter: N3 - CT - H1
Could not find angle parameter: N* - CC - CT
Could not find angle parameter: N* - CT - C
Could not find angle parameter: O - CC - N*
Could not find angle parameter: CC - CC - CM
Could not find angle parameter: CC - N* - CT
Could not find angle parameter: CC - CM - CA
Could not find angle parameter: CC - CM - HC
Could not find angle parameter: CC - CC - N*
Could not find angle parameter: CC - CC - O
Could not find angle parameter: NB - CC - N*
Could not find angle parameter: NB - CC - CM
Could not find angle parameter: NB - CC - CC
Could not find angle parameter: CC - NB - CC
Could not find angle parameter: CC - N* - CC
Could not find angle parameter: CC - N* - CT
Could not find angle parameter: CC - CT - N3
Could not find angle parameter: CC - CT - H1
Building proper torsion parameters.
** No torsion terms for CT-N*-CC-CT
** No torsion terms for N*-CC-CC-CM
** No torsion terms for O-CC-CC-CM
** No torsion terms for O-CC-N*-CT
** No torsion terms for CC-CC-CM-CA
** No torsion terms for CC-CC-CM-HC
** No torsion terms for CC-N*-CC-CT
** No torsion terms for CC-CC-N*-CT
** No torsion terms for NB-CC-N*-CC
** No torsion terms for NB-CC-N*-CT
** No torsion terms for NB-CC-CM-CA
** No torsion terms for NB-CC-CM-HC
** No torsion terms for NB-CC-CC-N*
** No torsion terms for NB-CC-CC-O
** No torsion terms for CC-N*-CC-CC
** No torsion terms for CC-N*-CC-O
Building improper torsion parameters.
total 4 improper torsions applied
Building H-Bond parameters.
Parameter file was not saved.
</pre>
This must then be set up in the Gaussian input file, two lines after the connectivity, in the following style:
=====Bonds=====
{| class="wikitable" border="1"
|-
!'''HrmStr1''': Harmonic stretch I (Amber 1): <math>ForceC(R-R_{eq})^{2}</math>
|-
|HrmStr1 Atom-type1 Atom-type2 ''ForceC'' <math>R_{eq}</math>

''ForceC'' Force constant

<math>R_{eq}</math> Equilibrium bond length
|}

=====Angles=====
{| class="wikitable" border="1"
|-
!'''HrmBnd1''': Harmonic bend (Amber 1): <math>ForceC(\theta-\theta_{eq})^{2}</math>
|-
|HrmBnd1 Atom-type1 Atom-type2 Atom-type3 ''ForceC'' <math>\theta_{eq}</math>

''ForceC'' Force constant (<math>Kcal mol^{-1}rad^{-2}</math>)

<math>\theta_{eq}</math> Equilibrium angle
|}

=====Torsions=====
{| class="wikitable" border="1"
|-
!'''AmbTrs''': Amber torsion (Amber 1): <math>\sum_{i=1}^{4} \frac{Mag_{i}[1+\cos(i\theta - POI(i+4))]}{N_{Paths}}</math>
|-
|AmbTrs Atom-type1 A-type2 A-type3 A-type4 ''PO1 PO2 PO3 PO4'' <math>Mag_{1}</math> <math>Mag_{2}</math> <math>Mag_{3}</math> <math>Mag_{4}</math> <math>N_{Paths}</math>

''PO1–PO4'' Phase offsets for <math>theta</math>: these may be set to 0 or 180: in the former case, they have no effect, in the latter, they have the sole effect of switching the sign of the '+1' coefficient in front of cos.

<math>Mag_{1}</math>-<math>Mag_{4}</math> <math>frac{V}{2}</math> magnitudes

<math>N_{Paths}</math> Number of paths. When zero or less, determined on-the-fly.
|}

Thus for the above example we would obtain a list that looks like:
<pre>
HrmStr1 CM HC
HrmStr1 CC N*
HrmStr1 CC O
HrmStr1 CC CM
HrmStr1 CC CC
HrmBnd1 CM C OH
HrmBnd1 CA C CM
HrmBnd1 HC CM CA
HrmBnd1 CM CA CM
HrmBnd1 CM CA CA
HrmBnd1 H1 CT N3
HrmBnd1 N* CC CT
HrmBnd1 N* CT C
HrmBnd1 O CC N*
HrmBnd1 CC CC CM
HrmBnd1 CC CM CA
HrmBnd1 CC CT N
HrmBnd1 CC CM HC
HrmBnd1 CC CC N*
HrmBnd1 CC CC O
HrmBnd1 NB CC N*
HrmBnd1 NB CC CM
HrmBnd1 NB CC CC
HrmBnd1 CC NB CC
HrmBnd1 CC N* CC
HrmBnd1 CC N* CT
HrmBnd1 CC CT N3
HrmBnd1 CC CT H1
AmbTrs CT N* CC CT
AmbTrs N* CC CC CM
AmbTrs O CC CC CM
AmbTrs O CC N* CT
AmbTrs CC CC CM CA
AmbTrs CC CC CM HC
AmbTrs CC N* CC CT
AmbTrs CC CC N* CT
AmbTrs NB CC N* CC
AmbTrs NB CC N* CT
AmbTrs NB CC CM CA
AmbTrs NB CC CM HC
AmbTrs NB CC CC N*
AmbTrs NB CC CC O
AmbTrs CC N* CC CC
AmbTrs CC N* CC O
</pre>
Now all that remains is to add the values to these parameters. To do this we go back to the files output when we ran the Redtools jobscript. Take the Mol-m1-o1-sm-mol2 file and open it in Gaussview. Change the PDB atom name and AMBER atom types of the .mol2 file in a text editor to those shown in the Gaussview atom list. This should be similar to the file below, obtained for GFP:

[[Image:Mol2_paramget.png|500 px|alt=File for getting AMBER parameters]]

Now we have this file we can obtain the missing parameters from the General AMBER Force Field using:
<pre>
parmchk -i Mol-m1-o1-sm.mol2 -f mol2 -o filename.frcmod
</pre>
where filename can be whatever you chose. The output of this file should now contain all the parameters required for the non-standard residues and the labels for the AMBER atom types should correspond directly to those output by ''saveamberparm''. For stretches and bends the numbers can be simply copied across, however, the torsions are a bit more complicated. An example of a torsion parameter from the .frcmod is shown below:
<pre>
H1-CT-C -O 1 0.800 0.000 -1.000 same as hc-c3-c -o
H1-CT-C -O 1 0.080 180.000 3.000 same as hc-c3-c -o
</pre>
Which must be put in the format

AmbTrs Atom-type1 A-type2 A-type3 A-type4 ''PO1 PO2 PO3 PO4'' <math>Mag_{1}</math> <math>Mag_{2}</math> <math>Mag_{3}</math> <math>Mag_{4}</math> <math>N_{Paths}</math>

Here H1, CT, C and O are the atom types; <math>N_{Paths}</math> is the second column; <math>mag_{i}</math> is the third column; ''POI'' is the fourth column; and the fifth column is the value of i/I. If there is a dash marker, this means that the next row is of the same torsion. The above example would translate then as:
<pre>
AmbTrs H1 CT C O 0 0 180 0 0.8000 0.0000 0.0800 0.0000 1.0
</pre>

We have now determined all the parameters for the AMBER calculation. In order to use them add ''amber=softfirst'' in the route section of the input file. A final point is that in the Gaussian input parameters, the atoms can be specified either way round (e.g. H1 CT C O or O C CT H1). These are equivalent and the input must be checked to ensure that each specification is unique, otherwise an error message will result. This happens even if the values are equal.

A example input file for GFP is shown here: [[Media:1W7S_01_SPE_amber.gjf|Gaussian AMBER input for GFP]]

{| class="wikitable" border="1"
|-
!Known Gaussian and Gaussview labelling problems
|-
|Hydroxyl protons are specified as HO in the AMBER atom types and in the parm96 force field parameters they have zero van der Waals radius. This results in Gaussian showing a warning that charged centres with zero van der Waals radii can collapse into a nearby oppositely charged centre, however, these centres should not have any radii associated with them.
|-
|Carbonyl oxygens are often specified as 'OM' in Gaussview, including any .com files it outputs. This is not recognized in any AMBER parameter sets, however, it is the same as 'O'. Any instances of 'OM' should be changed to 'O'. If this is not done the same problems of charged centres with no van der Waals radii as above occur.
|-
|The backbone nitrogens are often incorrectly labelled as 'N3', however, they should be labelled 'N'.
|-
|Aromatic carbon atoms are labelled 'CH', however, they should be 'CA'.
|}

==Constructing the ONIOM input==
We now have an Gaussian input file which will produce an AMBER calculation of the structure originally specified in the PDB database. This final section details how to progress from this point to an ONIOM input file. In order to check that everything is in order it may be worth running a single point AMBER calculation on the structure. This can be done using ''IOp(4/119=10)'' which will print out the force field parameters so that they can be checked. Providing there are no missing parameters this calculation should complete and an AMBER optimization can be carried out on this structure.

The converged AMBER structure can then be used to make the ONIOM file. A problem here is that there is no partial charge data upon opening the .log file in Gaussview, and the formatted checkpoint file loses all the PDB data. The easiest way to solve this then is to save the .log file as a .pdb file and then to open the .pdb file in Gaussview and save it as a .com file using Calculate→Gaussian Calculation Setup, although this will require re-entering the charges on the non-standard residue. Another option may be to use this [[Media:Onistruct.sh|script]] to add PDB data to the .com file obtained from the .fchk file. Although this is not necessary here, it is useful if you want to create an input with a geometry obtained from a calculation that used geom=check as the .log file loses all PDB data as well.

Whichever option you choose, open the file in Gaussview, specify the high level region using Edit→Edit Layer and then save as a .com file (using Calculate→Gaussian Calculation Setup as otherwise all MM charge data will be lost). The route section you use should look like this for mechanical embedding:
<pre>
#p opt oniom(b3lyp/6-31g(d):amber=softfirst) geom=connectivity
</pre>
or this for electronic embedding:
<pre>
#p opt oniom(b3lyp/6-31g(d):amber=softfirst)=embed geom=connectivity
</pre>

We now have a complete ONIOM input file: [[Media:1W7S_01_OPT_oniom_b3lyp_631gd_amber.gjf|1W7S_01_OPT_oniom_b3lyp_631gd_amber.gjf]]

If you try to run this calculation and get a missing parameter error, this is highly likely to be due to the fact that the parameters involving the link atoms are not present (these parameters could be obtained in the previous steps by using the actual model structure rather that the model without link atoms as done above and the above method should be modified to do this (perhaps using antechamber on a pdb file from Gaussview) in due course). If this happens, check that the missing parameters are not an indication of bad ONIOM partitioning (such as a link-atom replacing an electron-withdrawing group) and add the parameters by hand from the General AMBER Force Field parameter set (available from the AMBER website).

==Summary==
[[Image:Flow_chart_PDB_to_ONIOM.jpg|700 px|Flow-chart summarising procedure for turning PDB structure file to ONIOM input file.]]

Back to [https://wiki.ch.ic.ac.uk/wiki/index.php?title=ONIOM_for_biomolecules ONIOM for biomolecules]

Visualizing High Model Orbitals

2012-09-05T11:03:02Z

Lmt09: /* Aim */

==Visualizing High Model Orbitals==

===Aim===
When using ONIOM we are often interested in the effect of an environment on the model compound. For this reason it is often useful to visualize the orbitals of the high model region after a calculation. When trying to do this in Gaussview, however, it is found that the orbitals displayed are those of the low model. This is the reason why it was necessary to construct orbitals using ''guess=input'' in the ONIOM(CASSCF:AM1) examples. This tutorial explains how to access the orbitals after a calculation has been run, for example if we wish to localize orbitals to ensure the correct active space has been chosen in the previous example. In the case of QM:MM calculations it is possible to skip the punch orbitals step as there is only one set of orbitals on the checkpoint file; just get the model geometry and read the guess from the checkpoint file directly.

===System===
In this tutorial we examine the spurious transition state in the diels-alder cycloaddition between maleic anhydride and cyclohexadiene. We extract the high model orbitals and localize them to ensure that the correct active space has been chosen. This allows us to check that the choice of active space is not the cause of disagreement with the high real calculation, which indicates a symmetric transition state.

[[Image:TS_male_cyc.jpg|frame|Spurious transition state of maleic anhydride and cyclohexadiene]]

===Method===
====Punch Orbitals====
The first task is to obtain the orbitals if the high model in a format that can be read back in by Gaussian. This can be achieved using the ''punch=MO'' keyword but, in order to punch the high model orbitals we need to use a nonstandard route.

First we use ''testrt'' to obtain the standard route:
<pre>
-----------------------------------------------------------------
#p oniom(casscf(6,6)/sto-3g:hf/sto-3g) guess=read nosymm punch=MO
-----------------------------------------------------------------
1/38=1,52=2/1;
2/12=2,15=1,17=6,18=5,40=1/2;
1/38=1,52=2,53=3172/20;
3/6=3,11=9,16=1,25=1,30=1,116=-2/1,2,3;
4/5=1,17=6,18=6/1;
5/5=2,38=6/2;
6/7=2,8=2,9=2,10=2,28=1/1;
1/52=2,53=2032/20;
3/6=3,16=1,25=1,32=1,116=101/1,2,3;
4/5=1,17=6,18=6/1,5;
5/5=2,17=1000000,38=6/10;
6/7=2,8=2,9=2,10=2,28=1/1;
* 1/52=2,53=1022/20;
* 3/6=3,11=9,16=1,25=1,30=1,116=-2/1,2,3;
* 4/5=1,17=6,18=6/1;
* 5/5=2,38=6/2;
* 6/7=2,8=2,9=2,10=2,28=1/1;
* 1/52=2,53=3014/20;
99/5=1,9=1,10=32/99;
</pre>

The asterisked lines can be removed as they correspond to the low real system and adjust the last line to read <nowiki>99/10=32/99</nowiki>. We can now construct the input file to punch out the high-model orbitals. '''Remember''' to add ''cp fort.7 $WORK/$FLD/$FLNM.orbs'' to your jobscript file after the gaussian execution line as this file will contain the punched orbitals.

<pre>
%nprocshared=2
%mem=2000MB
%chk=/work/lmt09/PHD_Y2/MALA_CYHEX/ONIOM/macyhexdiene_S0_SPpunch_oniom_cas66_sto3g_hf_sto3g
# nonstd
1/38=1,52=2/1;
2/12=2,15=1,17=6,18=5,40=1/2;
1/38=1,52=2,53=3172/20;
3/6=3,11=9,16=1,25=1,30=1,116=-2/1,2,3;
4/5=1,17=6,18=6/1;
5/5=2,38=6/2;
6/7=2,8=2,9=2,10=2,28=1/1;
1/52=2,53=2032/20;
3/6=3,16=1,25=1,32=1,116=101/1,2,3;
4/5=1,17=6,18=6/1,5;
5/5=2,17=1000000,38=6/10;
6/7=2,8=2,9=2,10=2,28=1/1;
99/10=32/99;

#p oniom(casscf(6,6)/sto-3g:hf/sto-3g) guess=read nosymm punch=MO
Punch high model orbitals for localization

0 1 0 1 0 1
H 0 -0.26330500 -1.99941700 -1.21363800 H
C 0 0.78714900 1.41335300 0.46807900 H
C 0 1.26328700 1.66455000 -0.88162300 H
O 0 -1.92502500 -0.01205900 2.50907200 L
C 0 2.10729300 0.78800800 -1.50119200 H
H 0 0.93300900 2.55923500 -1.39147400 H
C 0 2.36471400 -0.53645200 0.68140600 L H 9 0.0000
C 0 -0.44009500 -1.08924500 0.83298900 H
C 0 1.08572800 0.18779800 1.13172100 H
C 0 -0.78359900 -1.36617800 -0.51617700 H
C 0 2.60670600 -0.49404400 -0.84831600 L H 5 0.0000
O 0 -2.42068300 0.13236600 0.27717100 L
C 0 -1.63821400 -0.28747100 1.36303200 L H 8 0.0000
H 0 2.12198000 -1.34178900 -1.33290700 L
O 0 -2.47630400 -0.31974200 -1.96661200 L
H 0 3.67299000 -0.60346100 -1.04369100 L
C 0 -1.93186300 -0.50956600 -0.89422500 L H 10 0.0000
H 0 0.08918800 2.11204200 0.91605600 H
H 0 -0.01783200 -1.85015200 1.47955400 H
H 0 0.96794500 0.21509600 2.21057000 H
H 0 2.35971100 -1.57164700 1.01714300 L
H 0 3.19920800 -0.04865300 1.18673800 L
H 0 2.42869300 0.97055600 -2.52117800 H

</pre>

====Obtain Geometry====
Now we have the orbitals for the high model in the .orbs file, however, if we wish to visualize these with gaussview we need to have them in a Gaussian output and so we need the geometry of the high model system. This can be done by taking the geometry output by the optimization and using the ''onlyinputfiles'' option.
<pre>
%nprocshared=1
%mem=800MB
%chk=/work/lmt09/PHD_Y2/MALA_CYHEX/ONIOM/macyhexdiene_S0_SPinput_oniom_cas66_sto3g_hf_sto3g
#p oniom(casscf(6,6)/sto-3g:hf/sto-3g)=onlyinputfiles nosymm

Input files

0 1 0 1 0 1
H 0 -0.26330500 -1.99941700 -1.21363800 H
C 0 0.78714900 1.41335300 0.46807900 H
C 0 1.26328700 1.66455000 -0.88162300 H
O 0 -1.92502500 -0.01205900 2.50907200 L
C 0 2.10729300 0.78800800 -1.50119200 H
H 0 0.93300900 2.55923500 -1.39147400 H
C 0 2.36471400 -0.53645200 0.68140600 L H 9 0.0000
C 0 -0.44009500 -1.08924500 0.83298900 H
C 0 1.08572800 0.18779800 1.13172100 H
C 0 -0.78359900 -1.36617800 -0.51617700 H
C 0 2.60670600 -0.49404400 -0.84831600 L H 5 0.0000
O 0 -2.42068300 0.13236600 0.27717100 L
C 0 -1.63821400 -0.28747100 1.36303200 L H 8 0.0000
H 0 2.12198000 -1.34178900 -1.33290700 L
O 0 -2.47630400 -0.31974200 -1.96661200 L
H 0 3.67299000 -0.60346100 -1.04369100 L
C 0 -1.93186300 -0.50956600 -0.89422500 L H 10 0.0000
H 0 0.08918800 2.11204200 0.91605600 H
H 0 -0.01783200 -1.85015200 1.47955400 H
H 0 0.96794500 0.21509600 2.21057000 H
H 0 2.35971100 -1.57164700 1.01714300 L
H 0 3.19920800 -0.04865300 1.18673800 L
H 0 2.42869300 0.97055600 -2.52117800 H

</pre>

====Localizing Orbitals====
The final part is to combine these two calculations to produce a Gaussian output with the orbitals on the high model. The first step is to copy and paste the high model input file from the ''onlyinputfiles'' output. '''Remember''' this is a single layer calculation so ensure that no ONIOM keywords are present. Once this is done add ''guess=cards'' to the route (note that in the example below the IOps have been removed as they are not necessary) and copy the .orbs file to the bottom of the input file. If localized orbitals are desired the vaarious keywords can be added here. This results in the following input:
[[Media:Male_loc.gjf]] 

We can now visualize the results which reveal a ''p-orbital'' on each carbon of the high model region, showing the correct active space has bee chosen.
[[Media:Male_loc.log]]

Visualizing High Model Orbitals

2012-09-05T11:02:08Z

Lmt09:

==Visualizing High Model Orbitals==

===Aim===
When using ONIOM we are often interested in the effect of an environment on the model compound. For this reason it is often useful to visualize the orbitals of the high model region after a calculation. When trying to do this in Gaussview, however, it is found that the orbitals displayed are those of the low model. This is the reason why it was necessary to construct orbitals using ''guess=input'' in the ONIOM(CASSCF:AM1) examples. This tutorial explains how to access the orbitals after a calculation has been run, for example if we wish to localize orbitals to ensure the correct active space has been chosen in the previous example. In the case of QM:MM calculations it is possible to skip the punch orbitals step as there is only one set of orbitals on the checkpoint file.

===System===
In this tutorial we examine the spurious transition state in the diels-alder cycloaddition between maleic anhydride and cyclohexadiene. We extract the high model orbitals and localize them to ensure that the correct active space has been chosen. This allows us to check that the choice of active space is not the cause of disagreement with the high real calculation, which indicates a symmetric transition state.

[[Image:TS_male_cyc.jpg|frame|Spurious transition state of maleic anhydride and cyclohexadiene]]

===Method===
====Punch Orbitals====
The first task is to obtain the orbitals if the high model in a format that can be read back in by Gaussian. This can be achieved using the ''punch=MO'' keyword but, in order to punch the high model orbitals we need to use a nonstandard route.

First we use ''testrt'' to obtain the standard route:
<pre>
-----------------------------------------------------------------
#p oniom(casscf(6,6)/sto-3g:hf/sto-3g) guess=read nosymm punch=MO
-----------------------------------------------------------------
1/38=1,52=2/1;
2/12=2,15=1,17=6,18=5,40=1/2;
1/38=1,52=2,53=3172/20;
3/6=3,11=9,16=1,25=1,30=1,116=-2/1,2,3;
4/5=1,17=6,18=6/1;
5/5=2,38=6/2;
6/7=2,8=2,9=2,10=2,28=1/1;
1/52=2,53=2032/20;
3/6=3,16=1,25=1,32=1,116=101/1,2,3;
4/5=1,17=6,18=6/1,5;
5/5=2,17=1000000,38=6/10;
6/7=2,8=2,9=2,10=2,28=1/1;
* 1/52=2,53=1022/20;
* 3/6=3,11=9,16=1,25=1,30=1,116=-2/1,2,3;
* 4/5=1,17=6,18=6/1;
* 5/5=2,38=6/2;
* 6/7=2,8=2,9=2,10=2,28=1/1;
* 1/52=2,53=3014/20;
99/5=1,9=1,10=32/99;
</pre>

The asterisked lines can be removed as they correspond to the low real system and adjust the last line to read <nowiki>99/10=32/99</nowiki>. We can now construct the input file to punch out the high-model orbitals. '''Remember''' to add ''cp fort.7 $WORK/$FLD/$FLNM.orbs'' to your jobscript file after the gaussian execution line as this file will contain the punched orbitals.

<pre>
%nprocshared=2
%mem=2000MB
%chk=/work/lmt09/PHD_Y2/MALA_CYHEX/ONIOM/macyhexdiene_S0_SPpunch_oniom_cas66_sto3g_hf_sto3g
# nonstd
1/38=1,52=2/1;
2/12=2,15=1,17=6,18=5,40=1/2;
1/38=1,52=2,53=3172/20;
3/6=3,11=9,16=1,25=1,30=1,116=-2/1,2,3;
4/5=1,17=6,18=6/1;
5/5=2,38=6/2;
6/7=2,8=2,9=2,10=2,28=1/1;
1/52=2,53=2032/20;
3/6=3,16=1,25=1,32=1,116=101/1,2,3;
4/5=1,17=6,18=6/1,5;
5/5=2,17=1000000,38=6/10;
6/7=2,8=2,9=2,10=2,28=1/1;
99/10=32/99;

#p oniom(casscf(6,6)/sto-3g:hf/sto-3g) guess=read nosymm punch=MO
Punch high model orbitals for localization

0 1 0 1 0 1
H 0 -0.26330500 -1.99941700 -1.21363800 H
C 0 0.78714900 1.41335300 0.46807900 H
C 0 1.26328700 1.66455000 -0.88162300 H
O 0 -1.92502500 -0.01205900 2.50907200 L
C 0 2.10729300 0.78800800 -1.50119200 H
H 0 0.93300900 2.55923500 -1.39147400 H
C 0 2.36471400 -0.53645200 0.68140600 L H 9 0.0000
C 0 -0.44009500 -1.08924500 0.83298900 H
C 0 1.08572800 0.18779800 1.13172100 H
C 0 -0.78359900 -1.36617800 -0.51617700 H
C 0 2.60670600 -0.49404400 -0.84831600 L H 5 0.0000
O 0 -2.42068300 0.13236600 0.27717100 L
C 0 -1.63821400 -0.28747100 1.36303200 L H 8 0.0000
H 0 2.12198000 -1.34178900 -1.33290700 L
O 0 -2.47630400 -0.31974200 -1.96661200 L
H 0 3.67299000 -0.60346100 -1.04369100 L
C 0 -1.93186300 -0.50956600 -0.89422500 L H 10 0.0000
H 0 0.08918800 2.11204200 0.91605600 H
H 0 -0.01783200 -1.85015200 1.47955400 H
H 0 0.96794500 0.21509600 2.21057000 H
H 0 2.35971100 -1.57164700 1.01714300 L
H 0 3.19920800 -0.04865300 1.18673800 L
H 0 2.42869300 0.97055600 -2.52117800 H

</pre>

====Obtain Geometry====
Now we have the orbitals for the high model in the .orbs file, however, if we wish to visualize these with gaussview we need to have them in a Gaussian output and so we need the geometry of the high model system. This can be done by taking the geometry output by the optimization and using the ''onlyinputfiles'' option.
<pre>
%nprocshared=1
%mem=800MB
%chk=/work/lmt09/PHD_Y2/MALA_CYHEX/ONIOM/macyhexdiene_S0_SPinput_oniom_cas66_sto3g_hf_sto3g
#p oniom(casscf(6,6)/sto-3g:hf/sto-3g)=onlyinputfiles nosymm

Input files

0 1 0 1 0 1
H 0 -0.26330500 -1.99941700 -1.21363800 H
C 0 0.78714900 1.41335300 0.46807900 H
C 0 1.26328700 1.66455000 -0.88162300 H
O 0 -1.92502500 -0.01205900 2.50907200 L
C 0 2.10729300 0.78800800 -1.50119200 H
H 0 0.93300900 2.55923500 -1.39147400 H
C 0 2.36471400 -0.53645200 0.68140600 L H 9 0.0000
C 0 -0.44009500 -1.08924500 0.83298900 H
C 0 1.08572800 0.18779800 1.13172100 H
C 0 -0.78359900 -1.36617800 -0.51617700 H
C 0 2.60670600 -0.49404400 -0.84831600 L H 5 0.0000
O 0 -2.42068300 0.13236600 0.27717100 L
C 0 -1.63821400 -0.28747100 1.36303200 L H 8 0.0000
H 0 2.12198000 -1.34178900 -1.33290700 L
O 0 -2.47630400 -0.31974200 -1.96661200 L
H 0 3.67299000 -0.60346100 -1.04369100 L
C 0 -1.93186300 -0.50956600 -0.89422500 L H 10 0.0000
H 0 0.08918800 2.11204200 0.91605600 H
H 0 -0.01783200 -1.85015200 1.47955400 H
H 0 0.96794500 0.21509600 2.21057000 H
H 0 2.35971100 -1.57164700 1.01714300 L
H 0 3.19920800 -0.04865300 1.18673800 L
H 0 2.42869300 0.97055600 -2.52117800 H

</pre>

====Localizing Orbitals====
The final part is to combine these two calculations to produce a Gaussian output with the orbitals on the high model. The first step is to copy and paste the high model input file from the ''onlyinputfiles'' output. '''Remember''' this is a single layer calculation so ensure that no ONIOM keywords are present. Once this is done add ''guess=cards'' to the route (note that in the example below the IOps have been removed as they are not necessary) and copy the .orbs file to the bottom of the input file. If localized orbitals are desired the vaarious keywords can be added here. This results in the following input:
[[Media:Male_loc.gjf]] 

We can now visualize the results which reveal a ''p-orbital'' on each carbon of the high model region, showing the correct active space has bee chosen.
[[Media:Male_loc.log]]

Resgrp:utilities

2012-09-05T10:54:19Z

Lmt09: /* gregap */

This is a share place for utilities, scripts, and Linux tricks (DMT).

== Bash shell ==

=== jump command (DMT) ===

The space quota for the $HOME directory is considerably smaller that $WORK. So, one should save input files in $HOME and save all large files (i.e. output, checkpoint files,...) in $WORK. Then, it is quite convenient to branch the folder directories identically in $HOME and $WORK, so they mirror each other.

The following command "jump" changes directory from $HOME to $WORK and vice-versa assuming that these are identically branched (use dirsync for this). Just add the following lines to your .bashrc profile.

<pre>
#! /bin/bash

jump () {
unset dname
unset jname
dname=$(pwd)
jname=$(echo $dname | cut -f 2 -d '/')
if [ "$jname" = "work" ] ; then
jname="${dname/work/home}"
echo "cd $jname"
else
jname="${dname/home/work}"
echo "cd $jname"
fi
cd $jname
}
</pre>

=== opt command ===
Allows you to view convergence criteria and step number from a Gaussian Log file.
<pre>
#!/bin/bash

egrep ' out of| SCF Don| Converged| NO | YES| exceeded' $1 | grep -v '\\\\'
</pre>

=== freq command ===
Prints out the first few frequencies, both in the 3N coordinates and the 3N-6 coordinates.
#!/bin/bash

egrep ' Low frequencies --| Frequencies --' $1 | grep -v '\\\\' | head -5
</pre>

=== dirsync command ===
Makes two directories, one in home, one in work. Useful to ensure your two directory trees are equivalent (use jump to go between them). if either directory tree is missing a portion you will get an error from the shell. you can also add commands to automatically generate files in a particular folder, here I create a file to add notes and the jobscript file I use in the home directory.

<pre>
#!/bin/bash

# creates directory in both home and work and copies gdvh13_118.sh to home dir

var1=$1;
name=$(pwd)
mkdir $name/${var1}

dname=$(pwd)
jname=$(echo $dname | cut -f 2 -d '/')
if [ "$jname" = "work" ] ; then
jname="${dname/work/home}"
mkdir $jname/${var1}
else
jname="${dname/home/work}"
mkdir $jname/${var1}
fi
path=$(echo $dname | cut -f 4- -d '/')

cp /home/lmt09/bin/gdvh13_118.sh /home/lmt09/$path/${var1}/gdvh13_118.sh
touch /home/lmt09/$path/${var1}/${var1}notes.txt
</pre>

=== cas command ===
View the MCSCF iterations to check convergence and also check any orbital rotations.

<pre>
#! /bin/awk -f

/ Step number/{
print;
}
/WARNING/{
line=$0;
getline;
line=line $0 ;
print line;
}
/ITN/{
print;
}
</pre>

== C shell ==

=== jump command ===

Same as the bash version but works for the c shell.
<pre>
#!/bin/csh

unset dname
unset jname
set dname="$PWD"
set jname=(`echo $dname | cut -f 2 -d '/'`)
switch ( $jname )
case work:
set jname=(`echo $dname | sed 's/work/home/'`)
echo "cd $jname"
breaksw
case home:
set jname=(`echo $dname | sed 's/home/work/'`)
echo "cd $jname"
breaksw
endsw
cd $jname
exit
</pre>

== Perl script==

=== gregap ===
graphically displays gaussian results from the log file using gnuplot (so you will need to be using X forwarding).

1. Download file

2. tar -xvzf gregap.tar.gz

3. ./install.sh (ensure this is executable -> chmod 755 install.sh)

[[Media:gregap.tar.gz | gregap.tar.gz]]

== Contacts ==

Here there is useful information and people to get in touch with:

Chief Service Technician (Mr. Peter Sulsh): p.sulsh@imperial.ac.uk‎

ICT Service Desk: service.desk@imperial.ac.uk‎

Allowances (Mrs. Althea Hartley-Forbes): a.hartley-forbes@imperial.ac.uk

HPC Systems Support Specialist (Mr. Matt Harvey): m.j.harvey@imperial.ac.uk

Printing posters: get in touch with ‎Dr. Ian R. Gould (i.gould@imperial.ac.uk). More information: http://www3.imperial.ac.uk/people/i.gould/poster%20printing1

File:Gregap.tar.gz

2012-09-05T10:52:09Z

Lmt09: uploaded a new version of "File:Gregap.tar.gz"

contains gregap files

Resgrp:utilities

2012-09-05T10:51:23Z

Lmt09:

This is a share place for utilities, scripts, and Linux tricks (DMT).

== Bash shell ==

=== jump command (DMT) ===

The space quota for the $HOME directory is considerably smaller that $WORK. So, one should save input files in $HOME and save all large files (i.e. output, checkpoint files,...) in $WORK. Then, it is quite convenient to branch the folder directories identically in $HOME and $WORK, so they mirror each other.

The following command "jump" changes directory from $HOME to $WORK and vice-versa assuming that these are identically branched (use dirsync for this). Just add the following lines to your .bashrc profile.

<pre>
#! /bin/bash

jump () {
unset dname
unset jname
dname=$(pwd)
jname=$(echo $dname | cut -f 2 -d '/')
if [ "$jname" = "work" ] ; then
jname="${dname/work/home}"
echo "cd $jname"
else
jname="${dname/home/work}"
echo "cd $jname"
fi
cd $jname
}
</pre>

=== opt command ===
Allows you to view convergence criteria and step number from a Gaussian Log file.
<pre>
#!/bin/bash

egrep ' out of| SCF Don| Converged| NO | YES| exceeded' $1 | grep -v '\\\\'
</pre>

=== freq command ===
Prints out the first few frequencies, both in the 3N coordinates and the 3N-6 coordinates.
#!/bin/bash

egrep ' Low frequencies --| Frequencies --' $1 | grep -v '\\\\' | head -5
</pre>

=== dirsync command ===
Makes two directories, one in home, one in work. Useful to ensure your two directory trees are equivalent (use jump to go between them). if either directory tree is missing a portion you will get an error from the shell. you can also add commands to automatically generate files in a particular folder, here I create a file to add notes and the jobscript file I use in the home directory.

<pre>
#!/bin/bash

# creates directory in both home and work and copies gdvh13_118.sh to home dir

var1=$1;
name=$(pwd)
mkdir $name/${var1}

dname=$(pwd)
jname=$(echo $dname | cut -f 2 -d '/')
if [ "$jname" = "work" ] ; then
jname="${dname/work/home}"
mkdir $jname/${var1}
else
jname="${dname/home/work}"
mkdir $jname/${var1}
fi
path=$(echo $dname | cut -f 4- -d '/')

cp /home/lmt09/bin/gdvh13_118.sh /home/lmt09/$path/${var1}/gdvh13_118.sh
touch /home/lmt09/$path/${var1}/${var1}notes.txt
</pre>

=== cas command ===
View the MCSCF iterations to check convergence and also check any orbital rotations.

<pre>
#! /bin/awk -f

/ Step number/{
print;
}
/WARNING/{
line=$0;
getline;
line=line $0 ;
print line;
}
/ITN/{
print;
}
</pre>

== C shell ==

=== jump command ===

Same as the bash version but works for the c shell.
<pre>
#!/bin/csh

unset dname
unset jname
set dname="$PWD"
set jname=(`echo $dname | cut -f 2 -d '/'`)
switch ( $jname )
case work:
set jname=(`echo $dname | sed 's/work/home/'`)
echo "cd $jname"
breaksw
case home:
set jname=(`echo $dname | sed 's/home/work/'`)
echo "cd $jname"
breaksw
endsw
cd $jname
exit
</pre>

== Perl script==

=== gregap ===
graphically displays gaussian results from the log file using gnuplot (so you will need to be using X forwarding).

1. Download file

2. tar -xvzf gregap.tar.gz

3. ./install.sh (ensure this is executable -> chmod 755 install.sh)

[[Media:gregap_1-0.tar.gz | gregap.tar.gz]]

== Contacts ==

Here there is useful information and people to get in touch with:

Chief Service Technician (Mr. Peter Sulsh): p.sulsh@imperial.ac.uk‎

ICT Service Desk: service.desk@imperial.ac.uk‎

Allowances (Mrs. Althea Hartley-Forbes): a.hartley-forbes@imperial.ac.uk

HPC Systems Support Specialist (Mr. Matt Harvey): m.j.harvey@imperial.ac.uk

Printing posters: get in touch with ‎Dr. Ian R. Gould (i.gould@imperial.ac.uk). More information: http://www3.imperial.ac.uk/people/i.gould/poster%20printing1

Resgrp:utilities

2012-05-18T10:10:16Z

Lmt09: /* GREGAP */

This is a share place for utilities, scripts, and Linux tricks (DMT).

== Bash shell ==

=== jump command (DMT) ===

The space quota for the $HOME directory is considerably smaller that $WORK. So, one should save input files $HOME and save all large files (i.e. output, checkpoint files,...) in $WORK. Then, it is quite convenient to branch the folder directories identically in $HOME and $WORK, so they mirror each other.

The following command "jump" changes directory from $HOME to $WORK and vice-versa assuming that these are identically branched (use dirsync for this). Just add the following lines to your .bashrc profile.

<pre>
#! /bin/bash

jump () {
unset dname
unset jname
dname=$(pwd)
jname=$(echo $dname | cut -f 2 -d '/')
if [ "$jname" = "work" ] ; then
jname="${dname/work/home}"
echo "cd $jname"
else
jname="${dname/home/work}"
echo "cd $jname"
fi
cd $jname
}
</pre>

=== opt command ===
Allows you to view convergence criteria and step number from a Gaussian Log file.
<pre>
#!/bin/bash

egrep ' out of| SCF Don| Converged| NO | YES| exceeded' $1 | grep -v '\\\\'
</pre>

=== freq command ===
Prints out the first few frequencies, both in the 3N coordinates and the 3N-6 coordinates.
#!/bin/bash

egrep ' Low frequencies --| Frequencies --' $1 | grep -v '\\\\' | head -5
</pre>

=== dirsync command ===
Makes two directories, one in home, one in work. Useful to ensure your two directory trees are equivalent (use jump to go between them). if either directory tree is missing a portion you will get an error from the shell. you can also add commands to automatically generate files in a particular folder, here I create a file to add notes and the jobscript file I use in the home directory.

<pre>
#!/bin/bash

# creates directory in both home and work and copies gdvh13_118.sh to home dir

var1=$1;
name=$(pwd)
mkdir $name/${var1}

dname=$(pwd)
jname=$(echo $dname | cut -f 2 -d '/')
if [ "$jname" = "work" ] ; then
jname="${dname/work/home}"
mkdir $jname/${var1}
else
jname="${dname/home/work}"
mkdir $jname/${var1}
fi
path=$(echo $dname | cut -f 4- -d '/')

cp /home/lmt09/bin/gdvh13_118.sh /home/lmt09/$path/${var1}/gdvh13_118.sh
touch /home/lmt09/$path/${var1}/${var1}notes.txt
</pre>

=== cas command ===
View the MCSCF iterations to check convergence and also check any orbital rotations.

<pre>
#! /bin/awk -f

/ Step number/{
print;
}
/WARNING/{
line=$0;
getline;
line=line $0 ;
print line;
}
/ITN/{
print;
}
</pre>

== C shell ==

=== jump command ===

Same as the bash version but works for the c shell.
<pre>
#!/bin/csh

unset dname
unset jname
set dname="$PWD"
set jname=(`echo $dname | cut -f 2 -d '/'`)
switch ( $jname )
case work:
set jname=(`echo $dname | sed 's/work/home/'`)
echo "cd $jname"
breaksw
case home:
set jname=(`echo $dname | sed 's/home/work/'`)
echo "cd $jname"
breaksw
endsw
cd $jname
exit
</pre>

== Perl script==

=== GREGAP ===
graphically displays gaussian results from the log file using gnuplot (so you will need to be using X forwarding). You will need to change some of the pathnames in the file.

[[Media:gregap_1-0.txt | gregap.pl]]

== Contacts ==

Here there is useful information and people to get in touch with:

Chief Service Technician (Mr. Peter Sulsh): p.sulsh@imperial.ac.uk‎

ICT Service Desk: service.desk@imperial.ac.uk‎

Allowances (Mrs. Althea Hartley-Forbes): a.hartley-forbes@imperial.ac.uk

HPC Systems Support Specialist (Mr. Matt Harvey): m.j.harvey@imperial.ac.uk

Printing posters: get in touch with ‎Dr. Ian R. Gould (i.gould@imperial.ac.uk). More information: http://www3.imperial.ac.uk/people/i.gould/poster%20printing1

2011-07-01T15:28:40Z

Lmt09:

File:Male loc.gjf

2011-07-01T15:28:28Z

Lmt09:

File:TS male cyc.jpg

2011-07-01T15:18:23Z

Lmt09:

Visualizing High Model Orbitals

2011-07-01T15:17:56Z

Lmt09: New page: ==Visualizing High Model Orbitals== ===Aim=== When using ONIOM we are often interested in the effect of an environment on the model compound. For this reason it is often useful to visuali...

==Visualizing High Model Orbitals==

===Aim===
When using ONIOM we are often interested in the effect of an environment on the model compound. For this reason it is often useful to visualize the orbitals of the high model region after a calculation. When trying to do this in Gaussview, however, it is found that the orbitals displayed are those of the low model. This is the reason why it was necessary to construct orbitals using ''guess=input'' in the ONIOM(CASSCF:AM1) examples. This tutorial explains how to access the orbitals after a calculation has been run, for example if we wish to localize orbitals to ensure the correct active space has been chosen in the previous example.

===System===
In this tutorial we examine the spurious transition state in the diels-alder cycloaddition between maleic anhydride and cyclohexadiene. We extract the high model orbitals and localize them to ensure that the correct active space has been chosen. This allows us to check that the choice of active space is not the cause of disagreement with the high real calculation, which indicates a symmetric transition state.

[[Image:TS_male_cyc.jpg|frame|Spurious transition state of maleic anhydride and cyclohexadiene]]

===Method===
====Punch Orbitals====
The first task is to obtain the orbitals if the high model in a format that can be read back in by Gaussian. This can be achieved using the ''punch=MO'' keyword but, in order to punch the high model orbitals we need to use a nonstandard route.

First we use ''testrt'' to obtain the standard route:
<pre>
-----------------------------------------------------------------
#p oniom(casscf(6,6)/sto-3g:hf/sto-3g) guess=read nosymm punch=MO
-----------------------------------------------------------------
1/38=1,52=2/1;
2/12=2,15=1,17=6,18=5,40=1/2;
1/38=1,52=2,53=3172/20;
3/6=3,11=9,16=1,25=1,30=1,116=-2/1,2,3;
4/5=1,17=6,18=6/1;
5/5=2,38=6/2;
6/7=2,8=2,9=2,10=2,28=1/1;
1/52=2,53=2032/20;
3/6=3,16=1,25=1,32=1,116=101/1,2,3;
4/5=1,17=6,18=6/1,5;
5/5=2,17=1000000,38=6/10;
6/7=2,8=2,9=2,10=2,28=1/1;
* 1/52=2,53=1022/20;
* 3/6=3,11=9,16=1,25=1,30=1,116=-2/1,2,3;
* 4/5=1,17=6,18=6/1;
* 5/5=2,38=6/2;
* 6/7=2,8=2,9=2,10=2,28=1/1;
* 1/52=2,53=3014/20;
99/5=1,9=1,10=32/99;
</pre>

The asterisked lines can be removed as they correspond to the low real system and adjust the last line to read <nowiki>99/10=32/99</nowiki>. We can now construct the input file to punch out the high-model orbitals. '''Remember''' to add ''cp fort.7 $WORK/$FLD/$FLNM.orbs'' to your jobscript file after the gaussian execution line as this file will contain the punched orbitals.

<pre>
%nprocshared=2
%mem=2000MB
%chk=/work/lmt09/PHD_Y2/MALA_CYHEX/ONIOM/macyhexdiene_S0_SPpunch_oniom_cas66_sto3g_hf_sto3g
# nonstd
1/38=1,52=2/1;
2/12=2,15=1,17=6,18=5,40=1/2;
1/38=1,52=2,53=3172/20;
3/6=3,11=9,16=1,25=1,30=1,116=-2/1,2,3;
4/5=1,17=6,18=6/1;
5/5=2,38=6/2;
6/7=2,8=2,9=2,10=2,28=1/1;
1/52=2,53=2032/20;
3/6=3,16=1,25=1,32=1,116=101/1,2,3;
4/5=1,17=6,18=6/1,5;
5/5=2,17=1000000,38=6/10;
6/7=2,8=2,9=2,10=2,28=1/1;
99/10=32/99;

#p oniom(casscf(6,6)/sto-3g:hf/sto-3g) guess=read nosymm punch=MO
Punch high model orbitals for localization

0 1 0 1 0 1
H 0 -0.26330500 -1.99941700 -1.21363800 H
C 0 0.78714900 1.41335300 0.46807900 H
C 0 1.26328700 1.66455000 -0.88162300 H
O 0 -1.92502500 -0.01205900 2.50907200 L
C 0 2.10729300 0.78800800 -1.50119200 H
H 0 0.93300900 2.55923500 -1.39147400 H
C 0 2.36471400 -0.53645200 0.68140600 L H 9 0.0000
C 0 -0.44009500 -1.08924500 0.83298900 H
C 0 1.08572800 0.18779800 1.13172100 H
C 0 -0.78359900 -1.36617800 -0.51617700 H
C 0 2.60670600 -0.49404400 -0.84831600 L H 5 0.0000
O 0 -2.42068300 0.13236600 0.27717100 L
C 0 -1.63821400 -0.28747100 1.36303200 L H 8 0.0000
H 0 2.12198000 -1.34178900 -1.33290700 L
O 0 -2.47630400 -0.31974200 -1.96661200 L
H 0 3.67299000 -0.60346100 -1.04369100 L
C 0 -1.93186300 -0.50956600 -0.89422500 L H 10 0.0000
H 0 0.08918800 2.11204200 0.91605600 H
H 0 -0.01783200 -1.85015200 1.47955400 H
H 0 0.96794500 0.21509600 2.21057000 H
H 0 2.35971100 -1.57164700 1.01714300 L
H 0 3.19920800 -0.04865300 1.18673800 L
H 0 2.42869300 0.97055600 -2.52117800 H

</pre>

====Obtain Geometry====
Now we have the orbitals for the high model in the .orbs file, however, if we wish to visualize these with gaussview we need to have them in a Gaussian output and so we need the geometry of the high model system. This can be done by taking the geometry output by the optimization and using the ''onlyinputfiles'' option.
<pre>
%nprocshared=1
%mem=800MB
%chk=/work/lmt09/PHD_Y2/MALA_CYHEX/ONIOM/macyhexdiene_S0_SPinput_oniom_cas66_sto3g_hf_sto3g
#p oniom(casscf(6,6)/sto-3g:hf/sto-3g)=onlyinputfiles nosymm

Input files

0 1 0 1 0 1
H 0 -0.26330500 -1.99941700 -1.21363800 H
C 0 0.78714900 1.41335300 0.46807900 H
C 0 1.26328700 1.66455000 -0.88162300 H
O 0 -1.92502500 -0.01205900 2.50907200 L
C 0 2.10729300 0.78800800 -1.50119200 H
H 0 0.93300900 2.55923500 -1.39147400 H
C 0 2.36471400 -0.53645200 0.68140600 L H 9 0.0000
C 0 -0.44009500 -1.08924500 0.83298900 H
C 0 1.08572800 0.18779800 1.13172100 H
C 0 -0.78359900 -1.36617800 -0.51617700 H
C 0 2.60670600 -0.49404400 -0.84831600 L H 5 0.0000
O 0 -2.42068300 0.13236600 0.27717100 L
C 0 -1.63821400 -0.28747100 1.36303200 L H 8 0.0000
H 0 2.12198000 -1.34178900 -1.33290700 L
O 0 -2.47630400 -0.31974200 -1.96661200 L
H 0 3.67299000 -0.60346100 -1.04369100 L
C 0 -1.93186300 -0.50956600 -0.89422500 L H 10 0.0000
H 0 0.08918800 2.11204200 0.91605600 H
H 0 -0.01783200 -1.85015200 1.47955400 H
H 0 0.96794500 0.21509600 2.21057000 H
H 0 2.35971100 -1.57164700 1.01714300 L
H 0 3.19920800 -0.04865300 1.18673800 L
H 0 2.42869300 0.97055600 -2.52117800 H

</pre>

====Localizing Orbitals====
The final part is to combine these two calculations to produce a Gaussian output with the orbitals on the high model. The first step is to copy and paste the high model input file from the ''onlyinputfiles'' output. '''Remember''' this is a single layer calculation so ensure that no ONIOM keywords are present. Once this is done add ''guess=cards'' to the route (note that in the example below the IOps have been removed as they are not necessary) and copy the .orbs file to the bottom of the input file. If localized orbitals are desired the vaarious keywords can be added here. This results in the following input:
[[Media:Male_loc.gjf]] 

We can now visualize the results which reveal a ''p-orbital'' on each carbon of the high model region, showing the correct active space has bee chosen.
[[Media:Male_loc.log]]