View source for Let's try with a smaller cycle

== Introduction ==

----
<p>
'''Advice break'''<br/>
First read the parts [https://www.ch.ic.ac.uk/wiki/index.php/Introduction Introduction], and [https://www.ch.ic.ac.uk/wiki/index.php/First_steps_with_ONIOM First steps with ONIOM]
</p>
----

<p>
Now try to run the same calculations with this new molecule :<br/>

[[Image:orbitals_flat.jpg|500px]] <br/>

(Benzene and a cycle of 8 carbones)<br/>
First, draw the molecule as above, then select the high and low layers as was done in the previous [http://www.ch.ic.ac.uk/wiki/index.php/First_steps_with_ONIOM tutorial]. The high layer should include the benzene ring, and the environment should be in the low layer. The cycle is smaller and so more stretched. It tends not to be flat. In order to get the right orbitals in the active space (without any mix between sigma and pi orbitals), it's better to force the benzene cycle to be flat during the orbitals section. Then, the geometry should be relaxed. 

</p>

----
<p>
'''Note Bene :''' <br/>
All these calculations were run with a development version of Gaussian (gdvh08_806).
</p>
----

== Orbitals ==

<p>
----
'''Advice break''' <br/>
Before beginning this part make sure that you read the part [https://www.ch.ic.ac.uk/wiki/index.php/ONIOM_for_excited_states ONIOM for excited states].
----
</p>

<p>
We need initial orbitals for the [http://www.gaussian.com/g_ur/k_casscf.htm CASSCF] computation. The active orbitals shall be chosen from these molecular orbitals.
</p>

<h4>Initial orbitals</h4>
<p>
In order to specify the correct orbitals for the CASSCF active space, they must be calculated them on the high model separately and then read these back into the ONIOM calculation. This is done using the ''onlyinputfiles'' option. It is important to use the ''nosymm'' keyword in the route section so that the molecular specification is consistent throughout.

<pre>
#p oniom(hf/sto-3g:am1)=onlyinputfiles pop=full nosymm
</pre>

</p>

<p>The output gives the three parts of the ONIOM calculation as three separate input files. The input file created by the ''inputfilesonly'' option, contains a number of [http://www.gaussian.com/g_tech/g_ur/m_molspec.htm ghost atoms]. It is not known the exact purpose of these, however, it appears that they act as 'place-holders' for the other atoms in the system that are outside the specific part of the model being calculated. Indeed, if they are not used then, when the orbitals are fed into the full ONIOM calculation, an error message is obtained.
</p>

<pre>
 #P Test IOp(2/15=1,5/32=2,5/38=1) RHF/STO-3G

 orbitals_flat
 Point  2 -- high level on model system.

     0     1
  C                                                -0.861573850000     -2.705360530000      0.001750920000
  C                                                 0.539826080000     -2.705430640000      0.002179930000
  C                                                 1.240587010000     -1.491817950000      0.001614740000
  C                                                 0.597351930000     -0.167438110000      0.003803130000
  C                                                -0.861451930000     -0.278065030000      0.000192820000
  C                                                -1.562212860000     -1.491677720000      0.000757740000
  H                                                 1.074779510000     -3.632104400000      0.002938500000
  H                                                 2.310586960000     -1.491871490000      0.001941870000
  H                                                -1.396405360000      0.648608730000     -0.000564350000
  H                                                -2.632212810000     -1.491624180000      0.000431570000
  H(Iso=12)                                        -1.295943091445     -3.457622827963      0.002100426539
  Bq-#1(Iso=1.00782504,Spin=1,GFac=2.792846)         -0.883254000000     -4.490611640000     -0.501545890000
  Bq-#1(Iso=1.00782504,Spin=1,GFac=2.792846)         -2.396418800000     -3.616885180000     -0.502505210000
  Bq-#1(Iso=12)                                    -1.719078100000     -4.188834530000      1.454108130000
  Bq-#1(Iso=1.00782504,Spin=1,GFac=2.792846)         -2.768405890000     -4.397214490000      1.473831030000
  Bq-#1(Iso=1.00782504,Spin=1,GFac=2.792846)         -1.183313270000     -5.015854910000      1.871113330000
  Bq-#1(Iso=12)                                    -1.120889690000     -3.151797470000      2.422779770000
  Bq-#1(Iso=1.00782504,Spin=1,GFac=2.792846)         -0.295044350000     -3.582882240000      2.949130780000
  Bq-#1(Iso=1.00782504,Spin=1,GFac=2.792846)         -1.757388370000     -2.640977520000      3.114756950000
  Bq-#1(Iso=12)                                    -0.491401960000     -2.062459720000      1.657126450000
  Bq-#1(Iso=1.00782504,Spin=1,GFac=2.792846)          0.419228440000     -2.410829260000      1.216341990000
  Bq-#1(Iso=1.00782504,Spin=1,GFac=2.792846)         -1.182658170000     -1.788160340000      0.887825010000
  H(Iso=12)                                         0.997112155446      0.603454641191      0.025971900794
  Bq-#1(Iso=1.00782504,Spin=1,GFac=2.792846)          0.548557320000      1.633329240000     -0.457719050000
  Bq-#1(Iso=1.00782504,Spin=1,GFac=2.792846)          2.046840460000      0.736917510000     -0.525984990000
  Bq-#1(Iso=12)                                     1.453825550000      1.276974390000      1.495605940000
  Bq-#1(Iso=1.00782504,Spin=1,GFac=2.792846)          2.481559650000      1.366224510000      1.779678200000
  Bq-#1(Iso=1.00782504,Spin=1,GFac=2.792846)          0.960159280000      2.214684190000      1.643571970000
  Bq-#1(Iso=12)                                     0.948223970000      0.163648320000      2.431806830000
  Bq-#1(Iso=1.00782504,Spin=1,GFac=2.792846)          1.573500390000     -0.585382130000      2.870991040000
  Bq-#1(Iso=1.00782504,Spin=1,GFac=2.792846)          0.444210660000      0.706048150000      3.204252960000
  Bq-#1(Iso=12)                                    -0.064162800000     -0.720106020000      1.679708890000
  Bq-#1(Iso=1.00782504,Spin=1,GFac=2.792846)         -0.757957550000     -0.453189660000      0.910095680000
  Bq-#1(Iso=1.00782504,Spin=1,GFac=2.792846)          0.810572600000     -1.149794240000      1.238000830000
</pre>

<p>
 The high model input can then be used to calculate initial orbitals at the STO-3G level. To do that, you can copy this section in a new input file. The first calculation is carried out on the STO-3G basis set using Hartree-Fock theory and then the size of the basis set is increased to 4-31G.

<h4>Visualise the orbitals and select the active space</h4>

The correct CASSCF active space must be specified in the relevant checkpoint file prior to the ONIOM calculation. This can be done using the orbitals from the RHF/4-31g single point energy calculation. By visualizing the orbitals in GaussView it can be determined that the relevant molecular orbitals are 18, 20, 21, 22, 23, 30. We can use the ''guess=alter'' keyword to swap orbitals to run the CASSCF/6-31g* single point energy calculation:

<pre>
#p cas(6,6)/6-31G(d) pop=full guess=(read,alter) geom=check nosymm

orbitals_cas_631gd

0 1

18 19
24 30
</pre>

{|
|-
| [[Image:active_space19.jpg|200px]]
| [[Image:active_space20.jpg|200px]]
| [[Image:active_space21.jpg|200px]]
|-
| [[Image:active_space22.jpg|200px]]
| [[Image:active_space23.jpg|200px]]
| [[Image:active_space24.jpg|200px]]
|}

<h4>Input and output files</h4>

[[Media:orbitals_flat.gjf]] <br/>
[[Media:orbitals_flat.log]] <br/>
[[Media:orbitals_flat_hf_sto3g.gjf]] <br/>
[[Media:orbitals_flat_hf_sto3g.log]] <br/>
[[Media:orbitals_flat_hf_431g.gjf]] <br/>
[[Media:orbitals_flat_hf_431g.log]] <br/>
[[Media:orbitals_cas_631gd.gjf]] <br/>
[[Media:orbitals_cas_631gd.log]] <br/>

----

'''Note Bene :''' <br/>
With the development version of Gaussian (gdvh08_806), the first output file seems to finish with an error. Nevertheless, the job is already done. So we can use the high model section to create the following input file all the same.

----

== Optimise the geometry of the ground state minimum ==

Now we can read the orbitals and run an optimization of the ground state minimum. As the benzene cycle doesn't need to be flat anymore, the geometry should be relaxed from now. This can be done by "cleaning" it with Gaussview.
</p>

<h4>Calculation</h4>

<p>
The energies of the different layers are calculated in the following order:<br/>
''Low Real''<br/>
''High Model'' <br/>
''Low Model'' <br/>
The low model and low real systems have not been calculated yet and so ''generate'' must be put so that these are calculated during the oniom calculation. For the high model system, the correctly ordered active space is now held in the checkpoint file of the CASSCF(6,6)/6-31G* single point energy calculation. This must now be read in using the ''guess=input'' keyword. We use an IOp(4/69=2) to make sure the orbitals are read from the checkpoint file and are not regenerated again.

<pre>
#p opt freq oniom(casscf(6,6)/6-31g(d):am1) nosymm guess=input geom=connectivity pop=full IOp(4/69=2)

''Molecular Specification''

generate

/work/mvacher/gdv/OniomCI/orbitals_cas_631gd.chk

generate
</pre>

Further information can be found at the bottom of the [http://www.gaussian.com/g_tech/g_ur/k_oniom.htm ONIOM user reference.] 
</p>

----

'''Nota Bene'''<br/>
With ''gdvh01'', the syntax is different. There is no extra blank lines.<br/>
<pre>
generate
/work/mvacher/gdv/OniomCI/orbitals_cas_631gd.chk
generate
</pre>
----

<h4>Results</h4>
<p>
The output should indicate that the high model initial orbitals have optimized within a few iterations, indicating that these have been read in correctly from the checkpoint file. The calculation should converge in 45 steps. The frequency calculation tell us the nature of this critical point. Check all frequencies are positive so this corresponds to the ground state minimum.<br/>
You can notice that you get the energy of the different "level of calculation" : 

<pre>
ONIOM: calculating energy.
ONIOM: gridpoint  1 method:  low   system:  model energy:     0.052031674589
ONIOM: gridpoint  2 method:  high  system:  model energy:  -230.759360919005
ONIOM: gridpoint  3 method:  low   system:  real  energy:    -0.011119897176
ONIOM: extrapolated energy =    -230.822512490770
</pre>

</p>

[[Image:s0_opt.jpg|500px]]

<h4>Input and output files</h4>
<p>
[[Media:s0_oniom_cas_631gd_am1.gjf]]<br/>
[[Media:s0_oniom_cas_631gd_am1.log]]


</p>

== Calculate the S1 Frank-Condon vertical excitation energy ==


<h4>Calculation</h4>
<p>
The vertical excitation energies can be calculated from this optimised geometry. The [http://www.gaussian.com/g_tech/g_ur/k_force.htm force] keyword can be used to provide information about the gradient of the potential energy surface at this geometry.<br/>

----

'''Keyword break'''<br/>
You must include '''nroot=x''' in the CAS keyword === > Calculations on excited states of molecular systems may be requested using the NRoot option. (Note that a value of 1 specifies the ground state, not the first excited state)<br/>

----

<pre>
 #p force oniom(casscf(6,6,nroot=2)/6-31g(d):am1) geom=check guess=read nosymm pop=full
</pre>
</p>

<h4>Results</h4>
<p>
You should get this energy in the output file. Ensure that the orbitals have converged within a few iterations so that we know they have been read in correctly. <br/>

<pre>
 ONIOM: calculating energy.
 ONIOM: gridpoint  1 method:  low   system:  model energy:     0.052034390477
 ONIOM: gridpoint  2 method:  high  system:  model energy:  -230.585348680238
 ONIOM: gridpoint  3 method:  low   system:  real  energy:    -0.011115630772
 ONIOM: extrapolated energy =    -230.648498701486
</pre>

The difference between the energy of the excited state and the ground state is:<br/>
ΔE = E<sup>S1</sup><sub>FC</sub> - E<sup>S0</sup><sub>min</sub> = 109.20 kcal/mol
</p>

<h4>Input and output fair</h4>
<p>
 [[Media:s0s1_FC_oniom_cas_631gd_am1.gjf]]<br/>
[[Media:s0s1_FC_oniom_cas_631gd_am1.log]]
</p>

== Optimise the geometry of the excited state minimum ==

<h4>Calculation</h4>
<p>
We will optimize the minimum on the S1 to be able to give the energy difference between the S1 minimum and the conical intersection.

<pre>
#p opt freq oniom(casscf(6,6,nroot=2)/6-31g(d):am1) geom=check guess=read nosymm pop=full
</pre>

</p>

<h4>Results</h4>
<p>
The geometry optimises in 12 steps with the following energies:

<pre>
ONIOM: calculating energy.
ONIOM: gridpoint  1 method:  low   system:  model energy:     0.069184790617
ONIOM: gridpoint  2 method:  high  system:  model energy:  -230.589044508969
ONIOM: gridpoint  3 method:  low   system:  real  energy:     0.001053820476
ONIOM: extrapolated energy =    -230.657175479110
</pre>

</p>

[[Image:s1_opt.jpg|500px]] <br/>

<h4>Input and output files</h4>
<p>
 [[Media:s1_oniom_cas_631gd_am1.gjf]]<br/>
[[Media:s1_oniom_cas_631gd_am1.log]]

</p>

== Find the S1/S0 conical intersections ==

----
<p>
'''Advice break''' <br/>
Before beginning this part make sure that you read the part [https://www.ch.ic.ac.uk/wiki/index.php/ONIOM_for_crossings ONIOM for crossings].
</p>
----

<p>
In order to find the conical intersection we will start from the Franck-Condon geometry. As was done in the [http://www.ch.ic.ac.uk/wiki/index.php/Resgrp:comp-photo-benzene-tutorial CASSCF tutorial], the planar structure of the benzene ring must be broken by moving one of the carbon atoms out of the plane as shown on the following picture. The molecule has C<sub>2</sub> symmetry so there are three possible carbon atoms to move out the plane. This suggests that there may be multiple different S1/S0 crossings that can be located by the conical intersection optimization. The conical intersection you'll find depend on which carbon atom you've moved out of the plane and towards which direction...   
</p>

[[Image:conical_guess.jpg|500px]]

<h4>Calculation</h4>

<p>
The fact that we are looking for a conical intersection is specified in the ''Opt=conical'' . 
<pre>
#p opt=conical oniom(casscf(6,6,nroot=2)/6-31g(d):am1) guess=read nosymm pop=full
</pre>
</p>

<h4>Results</h4>
<p>
The geometry optimizes in 11 steps. The optimized geometry is the following :
[[Image:conical.jpg|500px]]<br/>

And the energy of this conical intersection is: 

<pre>
ONIOM: calculating energy.
ONIOM: gridpoint  1 method:  low   system:  model energy:     0.218401034221
ONIOM: gridpoint  2 method:  high  system:  model energy:  -230.568625411086
ONIOM: gridpoint  3 method:  low   system:  real  energy:     0.147379795087
ONIOM: extrapolated energy =    -230.639646650220
</pre>

</p>

<p>
The difference between the energy of the conical intersection and S1 is:<br/>
ΔE = E<sup>S0/S1</sup><sub>CI</sub> - E<sup>S1</sup><sub>min</sub> = 11.00 kcal/mol
</p>

<h4>Compare with the benzene alone</h4>

<p>
Some similar calculations was realized on benzene (you will be able to find the results [https://www.ch.ic.ac.uk/wiki/index.php/Resgrp:comp-photo-benzene-tutorial here]), and on an other molecule ([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 Here])
</p>

<p>
So we can compare the difference <math>E_{crossing} - E_{S1}</math> between this 3 molecules. <br/>

What we hoped was that the fact that the cycle help to push the carbon out of the plan during the conical intersection optimization, so that the <math>E_{crossing} - E_{S1}</math> decrease compared to the benzene alone. <br/>

<h4>Input and output files</h4>

[[Media:s0s1_CI_oniom_cas_631gd_am1.gjf]]<br/>
[[Media:s0s1_CI_oniom_cas_631gd_am1.log]]


Back to [https://www.ch.ic.ac.uk/wiki/index.php/ONIOM_for_excited_states ONIOM for excited states]

Back to [https://www.ch.ic.ac.uk/wiki/index.php/ONIOM_for_crossings ONIOM for crossings]

Back to [https://www.ch.ic.ac.uk/wiki/index.php/ONIOM_tutorial_%28G03%29 ONIOM]