Labs 9 and 10: The use of solvation models and the ONIOM layered approach in Gaussian. In this lab we will consider two techniques that are very useful to model larger systems: the use of solvation models to mimic systems in solution, and the use of the ONIOM model to calculate different parts of the molecule at different levels of accuracy. Unfortunately they are mutually exclusive in the current implementation in Gaussian. Moreover the solvation model cannot be used in conjunction with semi-empirical methods. This is a limitation of Gaussian. There is no good reason why this would not be possible. A. Continuum solvation models in Gaussian03. The first type of calculation uses an implicit solvation model. It is known as the polarized continuum model (PCM) and was developed by Prof. Jacopo Tomasi from the University of Pisa. In this approach the interaction between a solvent and solute is modelled by a so- called Self-Consistent Reaction Field model. The molecule is placed inside a cavity, made up from the van der Waals sphere of the atoms in the molecule. Inside the cavity we have the permitivity of vacuum, outside the sphere we have a dielectricum with a dielectric constant that mimics the solvent of interest, for example water has 78.39 , while ether has 4.335 (in the Gaussian program). In the calculation the dipole of the solute interacts with the electric field of the solvent, and the electronic structure of the molecule is determined in the presence of this interaction. In the PCM model the electric field by the solvent is represented by apparent surface charges on the boundary of the cavity. This provides a very effective way of solving the equations. The use of a proper cavity is the most tricky part of the calculation, but it is handled automatically by Gaussian. Let us consider an example where solvation effects are important. We will consider various isomers of glycine, both in the gas phase and in water. Consider the following structures H O O O O - C O O C C H H C H H H N H C H N+ N H H H H H H H Isomer A Isomer B Zwitter Ion Let us first build these three versions of glycine in the gas phase, pre-optimize structures at the AM1 level, and then run full optimizations at the B3LYP / 6-31G(d) level. For more accurate structures we should really use a better basis set (e.g. 6-31++G(d,p)), but this takes too much time for our lab. Our optimization provides gas phase geometries. Next start from the B3LYP optimized gas phase structures and reoptimize the geometry, including effects of solvation. Under solvation select the default model (a variety of PCM) and select water as the solvent. Pull down the menu to see the other possibilities. The solvents are arranged in order of decreasing polarity (dielectric constant). Now we can optimize the three versions of glycine, mimicking the aqueous environment. The energies we get are not really good enough yet, because the basis set is lacking. For each of the optimized geometries and solvation models (gas phase or water) run single point energy calculations using the 6-31++G(d,p) basis set. This includes diffuse functions on all atoms (including H) and polarization functions on hydrogen also. This is in particularly important for the gas phase calculations. Moreover, on the keyword line put SCF=Tight. Gaussian prints out a warning, that this is needed with diffuse functions otherwise. Now we can fill in relative energies for these isomers in the gas phase and in solution. Assignments on solvation effects: 1. Finish the above exercise and calculate the relative energetics in KJ/mol for the above three structures both in vacuum and in water. Now repeat the set of calculations using ether as a solvent (both structure optimization at the B3LYP / 6-31G(d) level and higher accuracy single point calculations (B3LYP / 6-31++G(d,p) with SCF=Tight). Comment on your findings. 2. Also find the transition state for the hydrogen transfer between the structure of Isomer B and the Zwitter ion. I found it most convenient to first do a QST3 calculation at the AM1 level, calculating the force constants in every step. Using the optimized AM1 transition state you can then locate the transition state (using QST3) at the B3LYP/6- 31G(d) level, both in the gas phase and in solution (water solvent). Calculate the frequencies to check you have indeed located the transition state. Once you have found the transition state calculate the energetics at the higher level (B3LYP // 6-31++G(d,p), SCF=Tight). Now you can fill in the electronic energies for all of the structures. To make comparison with experiment we would need to also calculate vibrational frequencies to include the effects due to zeropoint motion and other thermal effects, but this time around let’s just say we did … B. ONIOM calculations in Gaussian03. Let us first do a warm-up problem to show the use of the ONIOM approach. Consider the calculation of the carbonyl stretch frequency in propanal, CH3CH2CHO. Low Level Link atom H H H C H H C H H C High level H O O In the (double layered) ONIOM approach we divide the molecule into two pieces. The region of most interest (containing the carbonyl group here), and the rest of the molecule. If we are interested in the CO bond length and frequency we can approximate the energy as E(ONIOM) = E(propanal / low level) + E(formaldehyde / high level) – E(formaldehyde / low level). We can use this compound formula for the energy to optimize the structure and to calculate the force constants. The positions of all the atoms are optimized for the complete structure, while the position of the link atom is somewhat arbitrary and uses a protocol developed in the group of Keji Morokuma, the inventor of ONIOM, at Emory University. Such a calculation is quite easy to set up in Gaussview. First draw propanal. Then select Multilayer Oniom model in the calculation set up under the method card. We now have to assign the atoms to the low and high level regions. This is done using the layer selection tool under the edit menu. You might first select all atoms in the low region (use set layer), and apply this selection. Next go on to the high layer, first initialize by select none and then select atoms by clicking on them and holding down the shift key. When you are done selecting the atoms use apply to make the changes permanent. In Gausview the atoms in different layers are rendered differently, so you can check your selection. Under view/display/format/molecule you can change the visualization. I prefer tubes in my low layer, in case of a double layer calculation. At this point we have selected the layers. In setting up the calculation we now have to specify the type of calculation in each layer. For our current application we may select the MP2/6-31G(d,p) for the high accuracy region, and HF/3-21G or AM1 for the low accuracy calculation. Gaussian will specify the link atom automatically for this simple calculation. It is useful to look at your input file (propanal.com) before running the calculation. You will see that each atom is assigned a level H or L. The carbon atom which is the link between the low and high accuracy regions has additional “H 1” on the input line. This specifies that in the high accuracy calculation the atom is replaced by H and it is bonded to atom 1 (the carbonyl carbon). Your numbering of the atoms may be different, it depends on your construction of the molecule. To see the atom numbering, use view/labels. This will show how the atoms are numbered in the input file. If the input file does not create the proper link atom, I think it is easiest to fix this in the input file directly. There are other possibilities to do this, which I think are more complicated. You can run the calculation and visualize frequencies as before. You will see that the carbonyl stretch does include some motion in your low accuracy region. We can increase the size of the calculation by selecting also the CH2 group to be part of the high accuracy region. The calculation will be more expensive, but you expect it to converge to the result we would get if we would do a complete MP2/6-31G(d,p) calculation. Let us finally extend the model one more step, and use a three-layer ONIOM model. We use formaldehyde in the high layer, using MP2 / 6-31G(d,p), then select the CH2 group for the second, medium layer, and will do a HF / 6-31G(d,p) calculation there, while finally we will use HF/3-21G calculations for the final CH3 group. This example is a little small for such a complicated approach of course, but it is easy to envisage how such an approach can be useful. It is definitely advantageous to make gradual changes in the computational level. Assignments on ONIOM: 1. The propanal molecule is small enough that we can also perform a full MP2 / 6- 31G(d,p) calculation. List the frequencies corresponding to the frequencies in the high accuracy region, and list the relevant geometrical parameters. Compare these quantities for all four (ONIOM) calculations that we did above, and identify the most appropriate ONIOM calculation. 2. Consider prototype electrophylic aromatic substitution reactions involving attack on benzene by either the Chlorine cation or the NO2 cation. These reactions involve intermediates (i.e. minima on the potential energy surface) of the following form where X = Cl, or X = NO2. We have seen these before. Don’t worry; no more transition states ... The intermediates have charge 1, and they are spin singlet states. a) Optimize the structures at the AM1 level and investigate the (Mulliken) charges. b) For both species substitute either the para or one of the meta hydrogens by a t-butyl group (one of the predefined groups in Gaussview). The t-butyl group is electron donating and you can expect that it will stabilize the intermediate if is attached to the carbon that has a positive charge in the parent cation under a). Let us investigate this by performing an ONIOM calculation. Select the methyl groups of the t-butyl group as your low level and the remainder (including the central carbon of the t-butyl group) as the high-level part of the ONIOM. We can use DFT/6-31G(d) for the high level and AM1 for the low level. Optimize your 4 structures (meta/para for either Cl or NO2 species), compare the energies of the para and meta substituted species and comment on the fact if this agrees with your expectations, based on the charge distributions seen under a). Also examine and compare the geometrical parameters involving the H-C-X moiety. c) For graduate students only: Carry out improved single point energy calculations at the ONIOM optimized geometries. We might use ONIOM with B3LYP/CC-PVTZ for the high accuracy region and 3-21G for the low accuracy region. Finally replace the t- butyl group by CH3 and carry out optimizations at the B3LYP/6-31G(d) level, followed by B3LYP/cc-PVTZ single point calculations (i.e. no ONIOM). Compare the geometries, and energetics between the methylated and t-butylated species. Also compare the charges on the most relevant atoms. I would not expect this substitution to have a huge impact, in particularly not on geometries. If the results differ substantially, it is more likely to be an error introduced by the ONIOM calculation on the t-butylated species. It is not entirely trivial to mix the proper computational levels and get accurate results from ONIOM. If you perform ONIOM calculations in actual research you should properly test things out, by making comparisons between gradually more accurate approaches.
Pages to are hidden for
"Labs 9 and 10_ The use of solvation models and the ONIOM layered "Please download to view full document