Resgrp:comp-photo-dyn/mctdh90.56dv/input
Creating Input for Dynamics Calculations
MCTDH (DD-vMCG) versions
- Development versions:
- mctdh90dv
- mctdh90.31dv
- mctdh90.56dv
The input required for the dynamics calculation is an optimised ground state structure and an optimised conical intersection in the appropriate orientation. [These two structures are used in the diabatic transformation. They need not necessarily correspond to the Franck-Condon and 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 optimised with state averaged orbitals. These orbitals must be "cleaned" and the high precision frequencies must have been calculated.
- Development version ( mctdh90.??dv): The only input required for the dynamics calculation is an optimised ground state structure and an optimised conical intersection. The rotation will be "automatic".
- 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
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
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.
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 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.
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. 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.)
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.) Use atom labels rather than numbers (otherwise the dd_generator will not work (see later)).
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.)
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.
Note
When using 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.
Keywords: CAS(4,4,nroot=2) STO-3G opt=(conical,tight) nosymm guess=read pop=full
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.