View source for Resgrp:comp-photo-dyn/mctdh90dev/input

= Creating Input for Dynamics Calculations =

== MCTDH versions ==

* Current version (mctdh90dev): The input required for the dynamics calculation is an optimised ground state structure and an optimised conical intersection in the appropriate orientation . The ground state structure must have been optimised with state averaged orbitals. These orbitals must be "cleaned" and the high precision frequencies must have been calculated.

* Development version (mctdh90.31dev): The only input required for the dynamics calculation is an optimised ground state structure and an optimised conical intersection (the rotation will be "automatic"). Again, the ground state structure must have been optimised with state averaged orbitals. These orbitals must be "cleaned" and finally high precision frequencies must be calculated.

* The main objective of this tutorial is learn how to run dynamics calculations, thus we will give you the optimised geometries in the appropiate rotation. (A description of how to rotate a structure will follow in due course)

== State Averaged Orbitals==

Discussion on the necessity of using state averaged orbitals to go here. NEED TO REPRESENT BOTH GROUND AND EXCITED STATES AT THE SAME TIME. I AM NOT CLEAR WHY YOU ARE USING STATE AVERAGING AT EQUILIBRIUM GEOMETRY. NORMALLY THIS WOULD NOT BE A SENSIBLE THING TO DO.

= Gaussian Calculations - CASSCF =

== '''1. Generate a starting geometry.''' ==

First of all, we need a set of initial orbitals for the CASSCF calculation, which can be obtained from a simple restricted Hartree-Fock (RHF/STO-3G) optimization. It is useful in systems with high symmetry to start with the molecule in the standard orientation and from a symmetrised geometry. This can be achieved by running a Gaussian calculation without the nosymm keyword and using the geometry returned by this calculation.

* NB You should use the structures given in this tutorial because they are  oriented appropriately.

[[Media:buta-s-orient.gjf]]

[[Media:buta-s-orient.log]]

== '''2. Identify the Orbitals Required for the CASSCF Calculation''' ==

Use the geometry from 1. Calculate the HF orbitals with a tight basis set e.g. STO-3G (so that it is easier to identify the relevant orbitals for CASSCF).

Keywords: '''HF/STO-3G pop=full nosymm geom=check'''

* Nosymm is required because every time a calculation is run on the molecule with a different symmetry Gaussian reorients the molecule for that particular point group and will destroy the active space. Nosymm disables this feature.
* Pop=full prints out the orbitals.
* Inspect the orbitals. Select those required for the CASSCF. 

[[Media:buta-hf.gjf]]

[[Media:buta-hf.log]]

[[Media:buta-hf.fchk]]

== '''3. Identify the State Averaged Orbitals for the Ground State''' ==

An orbital inspection using any visualization software before performing any CASSCF calculation is essential to identify which orbitals make up the chosen active space. In some cases it may be necessary to use '''guess=alter''' to change the order of the orbitals to make sure the optimisation (see step 4) proceeds with the correct orbitals. Gaussian assumes that the electrons specified in a CASSCF calculation come from the highest occupied orbitals in the initial guess and then takes the remaining orbitals as the lowest virtuals of this guess. 

* In this example, the relevant orbitals are the four pi orbitals, which contain four electrons (so the calculation is 4,4). The four orbitals that form the pi molecular orbitals are 14, 15, 16 and 17, thus we do not have to worry about swapping orbitals with the keword ''alter''.

Use the orbitals chosen in 2. (i.e. use the .chk from 2.) Use the final geometry from 2.

Keywords: '''CAS(4,4,nroot=2,stateaverage) STO-3G nosymm IOP(5/97=100,10/97=100) guess=read'''

Here, the weights of each state must be listed explicitly. The weights are 0.5 and 0.5. These must be specified by typing 0.5 followed by 7 spaces then 0.5 (on the same line) (The Gaussian program reads 10 characters (inc spaces) for the first weight and another 10 characters for the second weight. Thus the 7 spaces plus the three characters "0", "." and "5" make up the first 10 characters before the second weight is read.)

* The IOp (5/97=100,10/97=100) tells the program to switch the order of the two states. This is so that we calculate the ground state rather than the excited state. (Gaussian CASSCF calculations always perform optimisations on the highest state listed, so by switching the order we force optimisation of what is actually the lower state.)
* Nroot=2 specifies the ''first'' excited state so the orbitals will averaged over the ground state and the first excited state.

* Inspect the orbitals. Select those required for the CASSCF optimisation. The relevant orbitals are still 14, 15, 16 and 17.

[[Media:buta-orbs-sa.gjf]]

[[Media:buta-orbs-sa.log]]

[[Media:buta-orbs-sa.fchk]]

== '''4. Perform a CASSCF Optimisation on the Ground State with the State Averaged Orbitals''' ==

Use the orbitals chosen in 3. Use the .chk from 3.

Keywords: '''opt=conical pop=full CAS(4,4,nroot=2,stateaverage) STO-3G nosymm IOP(10/10=700005) IOP(5/97=100,10/97=100) geom=checkpoint guess=read'''

N.B. This is not a conical intersection optimisation! However, the keyword opt=conical specifies automatically that the populations of the two states will each be 0.5. Thus there is no need to specify the populations manually (which can be tedious). The initial '7' in IOp (10/10=700005) specifies that the calculation should be an optimisation to a minimum (rather than a conical intersection). The final '5' ensures that the orbital rotation derivative contributions from the CP-MC-SCF equations are included. (This is necessary and is usually the default. However, in the current version of Gaussian, the combination of '''opt''' and '''stateaverage''' excludes this).
Stateaverage is the default for opt=conical.
Nroot=2 is the default for opt=conical.
Do not include the keyword '''freq''' in the optimisation. (This combination of IOps does not allow '''opt''' and '''freq''' to be keywords in the same job.)

[[Media:buta-cas.gjf]]

[[Media:buta-cas.log]]

== '''5. "Clean" the Orbitals''' ==

In this case, all the orbitals in the active space are pure pi orbitals. Due to small numerical inaccuracies, the final orbitals may include small contributions from other (non-pi) orbitals. The aim of this step is to remove these contributions so that the active space consists (in this example) only of contributions from the C 2PY atomic orbitals. The four relevant orbitals are occupied orbitals 14 and 15 and virtual orbitals 16 and 17. In fact, they are already clean in this example. Verify this by inspection of the "Molecular Orbital Coefficients" in the population analysis. 

In other cases, multiple optimisations may be necessary to clean the orbitals:

Open the file in GaussView, and Edit/Symmetrize. Use this symmetrised geometry as a starting geometry for another optimsation as in 4 (including the keyword '''nosymm'''). 
Use the .chk from 4.
Re-run the optimisation exactly as in 4.

''Ideally here we include the keyword '''opt=(maxstep=1)'''. This should ensure that the optimisation takes the smallest possible step: we already know that the geometry is very close to the minimum and we need to prevent any significant re-orientation that might occur if large steps are taken. However, this does not seem to work (at least with g03).''

Repeat until the orbitals are clean.

== '''6. Calculate the (High Precision) Vibrational Modes for the Ground State ''with State Averaged Orbitals''''' ==

Use the .chk and the geometry from 5. The geometry must be specified explicitly. (With this combination of IOPs, Gaussian extracts the wrong geometry from the .chk file if '''geom=checkpoint''' is specified.)

Keywords: '''CAS(4,4,nroot=2) STO-3G freq=hpmodes guess=read IOP(5/17=41000200,10/10=700007) IOP(5/97=100,10/97=100)'''

The IOp (5/17=41000200,10/10=700007): includes the CP-MC-SCF correction (as above) (for a frequency calculation); and specifies state averaged orbitals. Do not include the keyword '''stateaverage''' with '''freq''' (the current version of Gaussian does not support this). Freq=hpmodes tells Gaussian to print the vibrational modes to 5 decimal places. As before, the IOp (5/97=100,10/97=100) forces calculation of the ground state. Do not use the keyword nosymm here. (Nosymm may prevent Gaussian from identifying the symmetries of the orbitals.)

[[Media:buta-hpfreq.gjf]]

[[Media:buta-hpfreq.log]]

== '''7. The Conical Intersection''' ==

One strategy for locating the conical intersection is to calculate the Frank-Condon point and run an optimisation from there, following (one of) the negative frequency(ies).

=== 7a. The Frank-Condon Point ===

Use the .chk from 4. Using the same orbitals perform a single point energy calculation for the first excited state.

Keywords: '''CAS(4,4,nroot=2,stateaverage) STO-3G pop=full nosymm guess=read geom=checkpoint'''

As in 3, the weights of each state must be listed explicitly. The weights are 0.5 and 0.5. These must be specified by typing 0.5 followed by 7 spaces then 0.5 (on the same line) ''(The Gaussian program reads 10 characters (inc spaces) for the first weight and another 10 characters for the second weight. Thus the 7 spaces plus the three characters "0", "." and "5" make up the first 10 characters before the second weight is read.)''

[[Media:buta-fc.gjf]]

[[Media:buta-fc.log]]

=== 7b. Optimise the Conical Intersection ===

Use the .chk from 7a. Distort the geometry to break the symmetry. 
A suggested distortion is to increase the C-C-C-C dihedral angle to 45° by twisting one of the terminal carbons out of the plane of the other atoms (be careful to move only the CH2 group, not the hydrogen on the adjacent carbon). Then increase the H-C-C-C dihedral angles of both the terminal CH2 groups to 70°, in a disrotatory direction. A conical intersection optimisation should now locate a conical intersection.

Keywords: '''CAS(4,4,nroot=2) STO-3G opt=conical nosymm guess=read pop=full'''

Stateaverage is the default for opt=conical.
Nroot=2 is the default for opt=conical.

[[Media:buta-dis14570ci.gjf]]

[[Media:buta-dis14570ci.log]]