Hands-on 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.

  1. To connect, the online command is: ssh -X login@styx.cbp.ens-lyon.fr
  2. put these lines in your .bashrc file:

    export g09root=/opt/64
    export GAUSS_SCRDIR=/tmp
    . \$g09root/g09/bsd/g09.profile

  3. connect to a node n (1...50) : ssh -X v22n
  4. 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 H2, the different steps needed to compute the Hartree-Fock wave function and look for the information we can extract from it. This will shed light on the limitation of an Hartree-Fock calculation.


1. First contact with Gaussian: Input and Output files

We shall first go back to basics and illustrate a few concepts on H2 (and get familiar with Gaussian).

We will study H2 with a minimal basis set (STO-3G), and illustrate RHF vs. UHF. The first thing is to create an input file H2.com containing the following lines:

#P HF STO-3G pop=regular test

Title: H2

0 1
H
H 1 0.7

  1. Submit the job and open the result file H2.log
  2. 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 ?

  3. From these links you will be able to find various informations. In particular, find:
    1. the number of MO, their symmetry
    2. the Hartree-Fock energy
    3. orbitals energies
    4. the kinetic energy
    5. the coefficients of the AO
  4. Can you write the HF wave function ?
  5. What is the value of J11 ?



2. Dissociating H2

To assess the strength of the H-H covalent bond, we need a reference for the system at infinite distance. It is composed of two not interacting hydrogen atoms.
  1. Get the energy of one hydrogen atom at the HF/STO-3G level of theory by readapting the previous input given for H2.
  2. Which method is used for the calculation, RHF of UHF ? Can you justify ?
  3. We now perform a so-called scan to retrieve the evolution of energy as a function of the inter-hydrogen distance d. We give the input H2scan.com that needs to be written in a scan.

    #P HF STO-3G 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 H-H distance by 0.2 A) each time from a starting distance of 0.5 A of the H-H 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.

  4. What do you observe for the energy at large distance ?
  5. We shall correct this unphysical behavior by following this procedure :
    • Please create two Gaussian input files with two Hydrogen atoms distant by 10 A. In one file specify a RHF calculation in the other one a UHF calculation, and run the calculations. What do you observe for the RHF calculation ? and for the UHF calculation ?
    • you can specify G09 to break the α, β and spatial symmetries: add guess=mix

      In this case the energy appears to be consistent as being the sum of the one of 2 isolated hydrogen atoms. However, what is the eignevalue of the S^2 operator ?

    • Try with a triplet state



3. An accurate geometry and λmax for C5H6NH2+ ? (DFT survey)

We now should consider the the tZt-penta-3,5-dieniminium cation that is a model for the rhodopsin. This tutorial is notably based on the two following publications publi-1 et publi-2. 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. some_text some_text

On pourra regarder ce lien pour une autre representation.

  1. We need a Z-matrix for this Retinal Minimal Model (2-cis-α-methyl-penta-2,4-dieniminium, 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 ?
  2. Optimize the cis geometry at the HF/6-31G* level of theory and then consider the optimizations using the BLYP, B3LYP, PBE (called PBE1PBE in G09), M06-2X, LC-BLYP, CAM-B3LYP 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 carbon-carbon or carbon-nitrogen 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 self-interaction error (see this paper for an in-depth 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 6-31+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" td-bench.log to plot λmax vs. the level of theory.
  • The EOM-CCSD 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 CAM-B3LYP prediction ? One can perform a T1 diagnostic (input, output) that suggests the need for multireference methods.



4. Multireference calculations : choosing an active space.

  • Re-optimize the cis-C5H6NH2+ cation at the HF/6-31G* level. Do not forget to save your checkpoint file by cp c5h6nh2.chk back.chk

    %mem=2000Mb
    %chk=c5h6nh2.chk
    #p hf/6-31G* 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) 6-31G* 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) 6-31G* 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 N2
We will compute the excitation energy of the N2 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 Single-Point calculations for the singlet ground state at the B3LYP level with the 6-31G* basis set, including the keywords gfinput iop(6/7=3).
  • Write down the molecular orbital energetic diagram of the N2 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 cc-pVTZ basis-set.
  • 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 cc-pVTZ basis-set, 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.
  1. To explore the PES, we need to decifer natural coordinate(s). Here it is legitimate to consider the two dihedral angles represented below.
    some_text
    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/6-31G* 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.

  2. Extract the (eight) geometries obtained after these two scans.
  3. 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".
  4. Before to refine our description, it is wise to get a first rude insight by inspecting the excited states with a "cheap" post Hartree-Fock method. Which is the cheapest ? The input will be as followed:

    #p METHOD 6-31G* 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.

  5. 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.
  6. If time allows, redo the study with a CASSCF(2,2) framework. Run these calculations and build up the variation of the S0 and S1 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