Resgrp:comp-photo-dyn/mctdh90.56dv Fulvene Tutorial/input
Creating Input for Dynamics Calculations
MCTDH (DD-vMCG) versions
- Development versions:
- mctdh90dv
- mctdh90.31dv
- mctdh90.56dv
The inputs required for the dynamic calculation are an optimised ground state structure and an optimized conical intersection structure in the appropriate orientation. [These two structures are used in the diabatic transformation. They do not necessarily need to correspond to the Franck-Condon and to the minimum conical intersection geometries. However, these structures are usually good places to start and so most of this tutorial will assume that these are the references that have been chosen]. The ground state structure must have been optimized with state averaged orbitals. These orbitals must be "cleaned" and the high precision frequencies must have been calculated.
- Development version ( mctdh90.??dv): The only inputs required for the dynamics calculation are an optimized ground state structure and an optimized conical intersection. The rotation will be "automatic".
- The main objective of this tutorial is to learn how to run dynamics calculations, thus we will give you the optimised geometries in the appropriate rotation. (A description of how to rotate a structure will follow in due course)
State Averaged Orbitals
State averaged orbitals are necessary in order to represent both ground and excited states at the same time in the space close to the conical intersection. The State Average algorithm allows the convergence of CASSCF when two or more states are close. A solution is found between the average of the involved orbitals if they were calculated independently for each state. At the Frank-Condon region, away from the intersection, this is not needed, but the code maintains this approximation over all the PES to avoid discontinuities.
Gaussian Calculations - CASSCF
Note: to use the dd_rotation method (see details for this method in the next section Rotation of the Conical Intersection) it is necessary to run all the calculations in this section with the gaussian development version gdvh01.
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) optimisation. 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.
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.
3. Identify the State Averaged Orbitals for the Ground State
An orbital inspection using any visualization software (e.g. GaussView) 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 the guess=alter keyword 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 six pi orbitals, which contain six electrons (so the calculation is 6,6). The six orbitals that form the pi molecular orbitals are 18, 20, 21,22, 23 and 24. Only orbital 18 is not included in the active space. Thus we need to swap orbitals 18 and 19 to get the correct active space. Use the keyword guess=alter .
Use the orbitals chosen in 2. (i.e. use the .chk from 2.) Use the final geometry from 2.
Keywords: CAS(6,6,nroot=2,stateaverage)/STO-3G nosymm IOP(5/97=100,10/97=100) guess=(read,alter)
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 be averaged over the ground state and the first excited state.
- Inspect the orbitals again after the calculation to check if the active space is correct.
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: CAS(6,6,nroot=2)/STO-3G nosymm opt=conical pop=full IOP(10/10=700005) IOP(5/97=100,10/97=100) geom=check 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.
- 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 some 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.)
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 six relevant orbitals are occupied orbitals 19, 20 and 21, and virtual orbitals 22, 23 and 24. In this example they are already clean. 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 4 (from 5 if you reoptimised the structure in order to clean the orbitals). 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(6,6,nroot=2)/STO-3G freq=hpmodes nosymm pop=full guess=read IOP(5/17=41000200,10/10=700007) IOP(5/97=100,10/97=100)
- 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.
- Use atom labels rather than numbers (otherwise the dd_generator will not work (see later)).
7. Optimise the Conical Intersection
Now we have to locate the conical intersection. To do this, we can run an optimisation from the Frank-Condon point (orbitals and geometry are the same we obtained with the CASSCF calculation in 4.) following (one of) the negative frequency(ies). Use the .chk from 4. A conical intersection optimisation should now locate a conical intersection.
- Stateaverage is the default for opt=conical.
- Nroot=2 is the default for opt=conical.
- To use the the dd_rotation method described in the next section (Rotation of the Conical Intersection) it may be convenient to use tighter condition to the optimization of the conical intersection at this stage (i.e. to use the keyword opt=tight).
The rational of this hack is that after rotation a second conical intersection optimization needs to be perfomed to obtain orbitals and branching space coordinates for the rotated geometry. These quantities cannot be obtained from a single point calculation. In the second optimization the less strict optimization conditions are used such that no actual geometry change is performed.
Keywords: CAS(6,6,nroot=2)/STO-3G opt=(conical,tight) nosymm guess=read pop=full