Handson 1st week RCTF. Illustrating HF, DFT and post HF methods in a nutshell.
Hands on : first week of RCTF Label Chimie Theorique.
Date: 29 and 30 October 2012,
P. Mignon and E. Dumont
This practical aims at using and applying the methods that were described during the lectures. We will use the Gaussian package (online help). Calculations will be launched on styx.
Visualization of geometries and inspection of molecular orbital will be done with molden.
 To connect, the online command is: ssh X login@styx.cbp.enslyon.fr
 put these lines in your .bashrc file:
export g09root=/opt/64
export GAUSS_SCRDIR=/tmp
. \$g09root/g09/bsd/g09.profile
 connect to a node n (1...50) : ssh X v22n
 you can perform a Gaussian calculation with: g09 < input.com > output.log
Here is a summary of what we will survey:
First, we will start with simple models such as H_{2}, the different steps needed to compute the HartreeFock wave function and look for the information we can extract from it. This will shed light on the limitation of an HartreeFock calculation.
1. First contact with Gaussian: Input and Output files
We shall first go back to basics and illustrate a few concepts on H_{2} (and get familiar with Gaussian).
We will study H_{2} with a minimal basis set (STO3G), and illustrate RHF vs. UHF.
The first thing is to create an input file H2.com containing the following lines:
#P HF STO3G pop=regular test
Title: H2
0 1
H
H 1 0.7
 Submit the job and open the result file H2.log
 Identify the various steps of the calculation with the reference to links (e.g.
/softs/gaussian/g03a02/g09/l1.exe, which you can search Link in the output file).
The specification of the links, available online will be useful to
understand the meaning of the link.
What are the sequential steps followed by G09 to compute the wave function and the energy of the system ?
 From these links you will be able to find various informations. In particular, find:
 the number of MO, their symmetry
 the HartreeFock energy
 orbitals energies
 the kinetic energy
 the coefficients of the AO
 Can you write the HF wave function ?
 What is the value of J_{11} ?
2. Dissociating H_{2}
To assess the strength of the HH covalent bond, we need a reference for the system at infinite distance. It is composed of two
not interacting hydrogen atoms.
 Get the energy of one hydrogen atom at the HF/STO3G level of theory by readapting the previous input given for H_{2}.
 Which method is used for the calculation, RHF of UHF ? Can you justify ?
 We now perform a socalled scan to retrieve the evolution of energy as a function of the interhydrogen distance d.
We give the input H2scan.com that needs to be written in a scan.
#P HF STO3G pop=regular test scan
Title: H2 scan
0 1
H
H 1 R1
R1 0.5 20 0.2
This will perform a scan (20 steps by increasing the HH distance by 0.2 A) each time from a starting distance of 0.5 A of the HH bond solving each time the RHF equations. Run the job and open the result file and look for Summary of the potential surface scan. A script is provided for more convenience.
 What do you observe for the energy at large distance ?
 We shall correct this unphysical behavior by following this procedure :
3. An accurate geometry and λ_{max} for C_{5}H_{6}NH_{2}^{+} ? (DFT survey)
We now should consider the the tZtpenta3,5dieniminium cation that is a model for the rhodopsin.
This tutorial is notably based on the two following publications publi1 et publi2. Let us first try to recover its vertical UV absorption.
To do so, we first need to have a good starting geometry (of course, one can obtain a good agreement with a given functional on a biaised geometry but this is certainly not a desirable feature).
The following publication adresses this problem in a more detailled way.
On pourra regarder ce
lien
pour une autre representation.
 We need a Zmatrix for this Retinal Minimal Model (2cisαmethylpenta2,4dieniminium, refered to as MPD hereafter).
It belongs to the family ofretinal protonated Schiff base (PSB) chromophores. Also consider the trans configurations : which one is the most stable ?
 Optimize the cis geometry at the HF/631G* level of theory and then consider the optimizations using the BLYP, B3LYP, PBE (called PBE1PBE in G09),
M062X, LCBLYP, CAMB3LYP and B2PLYP functionals and also the MP2 method.
 What are the requirements for the basis set here ?
 Why not using a dispersion correction for these functionals ?
 Classify these functionals along the Jacob's ladder.
 To obtain a reference geometry to judge the quality of the density functionals, MP2 may not be sufficient : why ?
Which other methods can you suggest ?
 Report the successive carboncarbon or carbonnitrogen distances. Which trend can you notice ?
This evolution is an example of the more general concept of bond length alternation (BLA), one of the useful metric for aromaticity. Double bond are shorter than single ones, but the difference of intercarbon distances is blurred because of the selfinteraction error
(see this paper for an indepth discussion).
We now move on to the calculation of a vertical transition for MPD. As excited state come into play, it is essential to include diffuse functions in the basis set : we will use 631+G**.
 We will be first using a TDDFT framework. The keyword is td. (The number of states generated depends on the size of the molecule of interest, three here.)
On the MP2 geometry, the calculations have been performed with the same subset of functionals than listed before. An input is given here, with
the corresponding output (the calculation takes some times...).
 Identify the MOs implied in the transition.
 Use the online command grep " nm" tdbench.log to plot λ_{max} vs. the level of theory.
 The EOMCCSD calculation provides a rather disappointing value (EOM stands for equation of motion). Suggest two reasons.
 What is still missing in our description ? Which calculation could provide us a legitimate reference value (cas1212.log, 4.352 eV) ? Is that for a good reason ?
 A legitimate question is to measure the error due to a bad geometry to the error caused by the electronic description in itself.
Can you say a word about that ?
 What would you say concerning the CAMB3LYP prediction ? One can perform a T1 diagnostic (input, output) that suggests the need for multireference methods.
4. Multireference calculations : choosing an active space.
 Reoptimize the cisC_{5}H_{6}NH_{2}^{+} cation at the HF/631G* level. Do not forget to save your checkpoint file by cp c5h6nh2.chk back.chk
%mem=2000Mb
%chk=c5h6nh2.chk
#p hf/631G* opt pop=full gfprint gfinput iop(1/6=60,1/7=600,5/6=6) nosymm
title
1 1
C
C 1 1.4
C 2 1.4 1 120.
C 3 1.4 2 120. 1 180.
C 4 1.4 3 120. 2 0.
N 5 1.4 4 120. 3 180.
H 1 1. 2 120. 3 0.
H 1 1. 2 120. 3 180.
H 2 1. 1 120. 8 0.
H 3 1. 2 120. 9 180.
H 4 1. 3 120. 10 0.
H 5 1. 4 120. 3 0.
H 6 1. 5 120. 12 0.
H 6 1. 5 120. 12 180.
We now can perform CASSCF calculation by increasing stepwisely the number of electrons and molecular orbitals of the active space. By default Gaussian will take the highest occupied and lowest unoccupied Molecular Orbitals as the active space. The associated keyword 'CAS(2,2)' is specified to perform a CASSCF calculation by taking into account 2 electrons and 2 moleuclar orbitals (1 occupied and 1 virtual).
 Run the calculation for the following input:
%mem=2000Mb
%chk=c5h6nh2.chk
#P CAS(2,2) 631G* pop=regular geom=allcheck guess=(read)
Start with 2 electrons allowed to be excited and 2 orbitals as the active space. Then try the different cases: (4 elec, 4 orb) and (6 elec, 6 orb).
In each case, check the CASSCF energy and determine which virtual orbitals are taken into account in the active space (see the orbital energy diagram from the optimization log file).
During a job the guess is changed, thus before running a job one has to take the initial checkpoint file. cp back.chk c5h6nh2.chk

Which orbitals are important and should be taken into account in the active space ? Which one would you include ?
 One can permute a given virtual orbital with another one more interesting or relevant. In the following example one allows N electrons and M orbitals to be included in
the active space and permute the virtual orbitals X and Y (Y is more interesting than X):
#P CAS(N,M) 631G* pop=regular geom=allcheck guess=(read,alter)
X, Y
After inspection of the orbital energy diagram, which active space and permutation would you propose to allow a better stabilization of the system ?
4bis. Excitation energy of N_{2}
We will compute the excitation energy of the N_{2} molecule from the singlet ground state to the triplet state B ^{3}Π_{g} which is the transition
σ_{g} > π_{g}.
Experimentally the value of 8.04 eV has been found. We will use different technics (crude and more sophisticated ones) to compute this excitations energy.
This exercice is inspired by this research article for those of you who want to go more deeply into details !
 Prepare the input files with a N...N distance of 1.0977 A (experimental value of the ground state) and launch SinglePoint calculations for the singlet ground state at the B3LYP level with the 631G* basis set, including the keywords gfinput iop(6/7=3).
 Write down the molecular orbital energetic diagram of the N_{2} ground state, including the 2s orbitals of the nitrogen atom (check the orbitals from molden and from the Gaussian output).
 Determine the desired transition from which we can determine B ^{3}Π_{g} ?
Which supplementary calculation should be performed to compute the corresponding excitation energy (check the electronic configuration in the supplementary calculation).
 Perform a CIS calculation to obtain excitations energies (ask to print the 10 first singlet and triplet states: link
with the ccpVTZ basisset.
 Check the output and look to the excitation energy corresponding to the good transition, (you will have to look to the symmetry of your RHF optimized orbitals to see the good numbers corresponding to the transtition: X > Y
 Perform CASSCF calculations with 10 electrons allowed to be in 10 orbitals with the ccpVTZ basisset, by using the same methodology as for B3LYP.
Conclusion: for all computed excited states energies, which method gives the best results ? Justify.
5. Looking for photoisomerization and searching a conical intersection...
We are now ready to focuss on the description of the photoisomerization of the MDP. We need some preliminary ground state calculations.
 To explore the PES, we need to decifer natural coordinate(s). Here it is legitimate to consider the two dihedral angles represented below.
We thus first need to generate four optimized geometries representative of the MDP cation torsions along φ_{1} and φ_{2}.
The header (for φ_{1}) is :
#p hf/631G* opt(addredu) pop=full gfprint gfinput iop(1/6=60,1/7=600,5/6=6) nosymm
1 1
[geom]
2 3 4 5 S 0. 3 30.
 Extract the (eight) geometries obtained after these two scans.
 Perform single point calculations in a dedicated directory to create a databank of checkpoint files. You will always need to inspect and note the orbital you "want".
 Before to refine our description, it is wise to get a first rude insight by inspecting the excited states with a "cheap" post HartreeFock method. Which is the cheapest ? The input will be as followed:
#p METHOD 631G* pop=full gfprint gfinput
iop(1/6=60,1/7=600,5/6=6) density=all geom(check) nosymm
From this job, you will be able to give : the absorption wavelength, the total energy of the system, the oscillator strength, the charge distribution, the dipole moment and the orbitals involved in the transition.
 We now perform the same CIS calculation along the angles φ_{1} and φ_{2} (0, 30, 60 and 90) and we will optimize the first excited state. You will need to add opt(addredu) in the keywords.
 If time allows, redo the study with a CASSCF(2,2) framework.
Run these calculations and build up the variation of the S_{0} and S_{1} states.
To find the mimimum energy. you will need to rotate the central dihedral by ~20 degrees.
geom00.com
geom030.com
geom060.com
geom090.com
geom300.com
geom600.com
geom900.com
Last update : Oct. 2012
