Page tree
Skip to end of metadata
Go to start of metadata

This part of the guide discusses how to:

  • Calculate the vertical excitation spectrum with TD-DFT
  • Carry out an excited state geometry optimization with TD-DFT 
  • Analyze the results from excited state calculations

We will use TURBOMOLE modules ridftescf, egrad, and the optimization script jobex. To follow this part of the guide, you need to the optimized structure of borazine from the previous part of the guide.

Table of contents

Setting up a vertical excitation calculation with escf

You should now be in the directory which has the TURBOMOLE inputs and results for borazine.

Creating a new subdirectory ex

It is usually more convenient to carry out the vertical excitation calculation in a subdirectory of the optimized ground state geometry. Create a new subdirectory by executing:

cpc ex

This copies the inputs necessary to run TURBOMOLE to a new subdirectory called ex (cpc = CoPy Control). Enter the new subdirectory and delete two files:

cd ex
rm hessapprox wherefrom

The hessapprox file contains the approximate Hessian from the previous geometry optimization. It would cause problems when we carry out an excited state optimization. The other file wherefrom is just an information file create by the cpc sript.

Setting up the escf calculation with define

Next, we add the vertical excitation setup into control with define. Execute

define

which will say

DATA WILL BE TAKEN FROM control BY DEFAULT

INPUT TITLE

No need to enter title, just press <Enter> 

Skip the geometry, basis set, and orbital menus by pressing <Enter> three times. Then, define will ask:

DO YOU WANT TO DELETE DATA GROUPS LIKE
$energy
$grad
$dipole
LEFT OVER FROM PREVIOUS CALCULATIONS ? DEFAULT(n)

We want to clean the ground state results from control, so answer y. Next, define asks three questions

delete data group $energy ? DEFAULT(y)
delete data group $grad ? DEFAULT(y)
delete data group $dipole ? DEFAULT(y)

Press <Enter> for each question and you will finally arrive in the general (method) menu. Execute

ex

to open the excited state calculation submenu:

Excited state submenu of define
 MAIN MENU FOR RESPONSE CALCULATIONS

 OPTION | STATUS | DESCRIPTION
 -------------------------------------------------------------------
 rpas   | off    | RPA SINGLET EXCITATIONS (TDHF OR TDDFT)
 ciss   | off    | TDA SINGLET EXCITATIONS (CI SINGLES)
 rpat   | off    | RPA TRIPLET EXCITATIONS (TDHF OR TDDFT)
 cist   | off    | TDA TRIPLET EXCITATIONS (CI SINGLES)
 polly  | off    | STATIC POLARIZABILITY
 dynpol | off    | DYNAMIC POLARIZABILITY
 single | off    | SINGLET STABILITY ANALYSIS
 triple | off    | TRIPLET STABILITY ANALYSIS
 nonrel | off    | NON-REAL STABILITY ANALYSIS

 ENTER <OPTION> TO SWITCH ON/OFF OPTION, * OR q TO QUIT

We want to calculate vertical singlet excitations with TDDFT. Execute

rpas

and then move ahead by pressing *. In the next step, we need to select the excited states that will be calculated:

Excited state selection in define
 STATE SELECTION MENU

 IRREP |  #STATES  | #SELECTED |
 ----------------------------------------------------------------
  a1'  |      230  |        0  |
  a2'  |      193  |        0  | Lz
  e'   |      423  |        0  | (x,y) (xx-yy,xy)
  a1"  |      100  |        0  |
  a2"  |      128  |        0  | z xz
  e"   |      228  |        0  | (Ly,Lx)

 SELECT IRREP AND NUMBER OF STATES
 ENTER ? FOR HELP, * OR Q TO QUIT, & TO GO BACK

The states are selected based on irreducible representations of the point group (D3h for borazine). Here a key point is that only irreps with the dipole operators (x,y,z) possess allowed electronic transitions. So, in the present case, let's select the 5 lowest energy vertical singlet excitations in the irreps e' and a2" by executing

e' 5
a2" 5

After this, you can close the state selection menu with *.

The next step is setting up the general options for the response calculations (escf):

General options for the response calculations submenu in define
 GENERAL OPTIONS FOR RESPONSE CALCULATIONS

 OPTION  | VALUE  | DESCRIPTION
 -----------------------------------------------------
 rpacor  |    200 | MEMORY (MB) (DEFAULT:   200)
 rpaconv |      5 | CONVERGENCE THRESHOLD (DEFAULT: 5)

 ENTER <OPTION> <VALUE> TO CHANGE OPTION, * OR q TO QUIT

Here it is important to increase the amount of memory for the excited state calculation by executing

rpacor 3000

The exact amount depends on your system resources, but for the computing clusters at Aalto CMAT, 3000 MB (per CPU core) is a good choice. Close the menu by executing *.

The next question from define is:

SET SCF DENSITY CONVERGENCE THRESHOLD $denconv TO 1d-7
(RECOMMENDED FOR RESPONSE CALCULATIONS) (y/n, DEFAULT:y)?

It is a good idea to calculate a bit more accurate wavefunction before the excited state calculation. Press <Enter> and define will note in the mos file that the MOs need to be converged again.

After this, define will again enter the general menu and you can finish by executing *.

Additional changes to the control file

Open the control file in an editor and you will see that the following keywords were added by define:

$scfinstab rpas
$soes
e' 5
a2" 5
$rpacor 3000
$denconv 1d-7

Increase also the SCF convergence criterion $scfconv from 6 to 7:

$scfconv   7

Add keyword $spectrum to get the vertical excitation spectrum in an convenient format (in nm units):

$spectrum nm

Other possible units are eV or rcm (for cm-1). Save the changes to the control file and we are ready to submit the job. 

Submitting a vertical excitation calculation

Before running the actual vertical excitation calculation with the escf module, we need to converge the MOs with the tighter convergence criteria ($denconv, $scfconv). Submit a single-point energy calculation with

jsub -mem 1G tm borazine_pbe0_svp_denconv ridft

Here we limit the memory reserved for the job to 1 GB. Otherwise jsub will interpret the $rpacor keyword in the control file and reserve too much.

After the job has finished, check that SCF converged OK. The file borazine_pbe0_svp_denconv.out should contain the line

convergence criteria satisfied after     3 iterations

Now, you can submit the actual escf calculation:

jsub -mem 4G tm borazine_pbe0_svp_ex escf

We ask for 4 GB of memory per CPU core because $rpacor was set to 3 GB and the job needs memory also in addition to $rpacor. Please note that for a small molecule such as borazine this is a ridiculously large amount of memory, but the purpose of the guide is to give as general guidelines as possible.

If you forgot to add the $spectrum keyword, jsub will add it when you submit the job.

Job submission commands for production jobs

When you are running larger production jobs in parallel, the submission commands on Aalto CMAT clusters are:

jsub -np 12 -smp mem 2G tm NNNN_denconv ridft

for the ridft single-point energy calculation and

jsub -np 12 -smp -mem 5G tm NNNN_ex escf

for the escf calculation. The -smp switch is important since it restricts the jobs on a single node (SMP parallel versions of ridft end escf will be used). 

Analyzing the results from a vertical excitation calculation

After escf has finished, it creates a file called spectrum. In this case, the file contains the excitation energies of the vertical singlet excitations and their oscillator strengths.

spectrum file for borazine
# Electronic excitation spectrum of , IRREP e'
# singlet excitations
#  excitation energy / nm  oscillator strength (length rep.)
 0.15493015277554E+03   0.72427546155199E+00
 0.12608130098241E+03   0.12821480233044E+00
 0.12486905622603E+03   0.70493035396661E-03
 0.12132361852177E+03   0.14599108976763E-02
 0.11614549361076E+03   0.18081358174180E-01
# Electronic excitation spectrum of , IRREP a2"
# singlet excitations
#  excitation energy / nm  oscillator strength (length rep.)
 0.15661992302118E+03   0.36809837579158E-01
 0.13417085320734E+03   0.45646785314809E-01
 0.11049090673867E+03   0.48308293702068E-03
 0.10967440403225E+03   0.13895172356267E-02
 0.10653787912137E+03   0.14187361804253E+00

The S0 → S1 excitation (a2") is at 157 nm (oscillator strength 0.037) and S0 → S2 (e') is at 155 nm (osc. strength 0.72). The larger the oscillator strength, the higher the probability of electronic absorption. In the case of borazine, experimental vacuum UV-Vis study shows an absorption maximum at 165 nm. The difference between 165 and 155 nm is big (0.5 eV), but we used a very small basis set. Increasing the basis set to def2-TZVP would already decrease the difference between experiment and theory to 0.28 eV.

Excited state density difference

Let's analyze the excited state by looking at electron density differences between the ground state and the excited state.

Edit the control file so that only the lowest energy excitation of the a2" irrep is included in $soes and $pointval keyword is added:

$soes
a2" 1
$pointval

Next, submit a excited state gradient calculation:

jsub -mem 4G tm borazine_pbe0_svp_s0_s1_ex egrad

The calculation will produce three files: coord.xyz (XYZ coordinates of the current geometry), td.plt (total electron density of the chosen excited state), and ed.plt (density difference between the ground state S0 and the chosen excited state S1.

Give ed.plt a more descriptive name:

mv ed.plt s0_s1_ex.plt

Calculate also the electron density difference plot for the S0 → S2 excitation with much higher oscillator strength. Change $soes keyword in control to

$soes
 e'      1

and submit the job:

jsub -mem 4G tm borazine_pbe0_svp_s0_s2_ex egrad

After the job has finished, rename ed.plt:

mv ed.plt s0_s2_ex.plt

Job submission commands for production jobs

When you are running larger production jobs in parallel, the excited state density difference submission command on Aalto CMAT clusters is:

jsub -np 12 -smp -mem 5G tm NNNN_s0_s1_ex egrad

Visualizing excited state density differences with VMD

For this task, you need to have VMD installed on your local workstation. In addition, copy this vmd.rc file. Replace the vmd.rc file located in the VMD installation folder or at least copy the parts after ### VMD-TM ### into your own vmd.rc file. 

Transfer the files coord.xyz, s0_s1_ex.plt and s0_s2_ex.plt to your local workstation.

Start VMD and run the following steps:

  1. File → New molecule
  2. Load coord.xyz
  3. Execute script akinit in the console to get better representation of the molecule (it will hide hydrogen atoms by default)
  4. File → New molecule
  5. Load the plt file with the density data
  6. Execute script akdens ex to illustrate the density. 

Blue color corresponds to volumes where the electron density increases during the S0 → Sn excitation.

Red color corresponds to volumes where the electron density decreases during the S0 → Sn excitation.

Setting up and submitting an excited state geometry optimization

Next, we will optimize the geometry of borazine for the lowest energy singlet (S1) and triplet (T1) excited states. Even though the S0 → S2 excitation had much larger oscillator strength than the S0 → S1 excitation, Kasha's rule states that normally emission (fluorescence or phosphorescence) occurs in appreciable yield only from the lowest excited state of a given multiplicity. Therefore, we will optimize the geometries of the Sand T1 excited states to study the emission energies.

The easiest way to start an excited state geometry optimization is to use the vertical excitation calculation as a starting point. Go to the directory ex and execute

cpc ../s1_opt

for the geometry optimization of the S1 state (the a2" 1 state we plotted above) and execute

cpc ../t1_opt 

for the geometry optimization of the T1 state (triplet version of the a2" 1 state).

Setting up the optimization of the S1 state

Go to the directory s1_opt. Now we need to tell TURBOMOLE that we want to optimize the geometry of the S1 excited state. Set the keyword $soes to:

$soes
a2" 1

and remove the $pointval keyword (creating the plot at every step of the geometry optimization would waste computational resources).

Note also that the $spectrum keyword should remain in the control file, this allows to easily monitor how the excitation energy develops during the optimization.

The excited state geometry optimization can be submitted with

jsub -mem 4G tm borazine_pbe0_s1_opt jobex -ri -ex

Note the -ex switch to jobex. This switch tells jobex to use egrad module to calculate the gradients.

Setting up the optimization of the T1 state

Go to the directory t1_opt. Now we need to tell TURBOMOLE that we want to optimize the geometry of the T1 excited state. Change the keyword $scfinstab to:

$scfinstab rpat

so that we optimize the triplet state instead of singlet.

Make sure that the keyword $soes is as follows:

$soes
a2" 1

and remove the $pointval keyword.

The excited state geometry optimization can be submitted with

jsub -mem 4G tm borazine_pbe0_t2_opt jobex -ri -ex

Job submission commands for production jobs

When you are running larger production jobs in parallel, the submission command for an excited state geometry optimization on Aalto CMAT clusters is:

jsub -np 12 -smp -mem 5G tm NNNN_s1_opt jobex -ri -ex

Analyzing the results from an excited state geometry optimization

Create a synthetic Gaussian log of the geometry optimization:

t2g -o > borazine_pbe0_svp_s1_opt.log

The final step is to prepare excited state density difference plots for the emission processes.

Excited state density difference for the emission

Both in the s1_opt and t1_opt directories, add $pointval to control file and submit one more egrad job :

jsub -mem 4G tm borazine_pbe0_svp_s1_s0_em egrad

for the S1 → S0 emission or

jsub -mem 4G tm borazine_pbe0_svp_t1_s0_em egrad

for the T1 → S0 emission.

After the calculations have finished, give ed.plt a more descriptive name in both directories:

mv ed.plt s1_s0_em.plt

and 

mv ed.plt t1_s0_em.plt

Job submission commands for production jobs

When you are running larger production jobs in parallel, the excited state density difference submission command on Aalto CMAT clusters is:

jsub -np 12 -smp -mem 5G tm NNNN_s1_s0_em egrad

for the S1 → S0 emission and

jsub -np 12 -smp -mem 5G tm NNNN_t1_s0_em egrad

for the T1 → S0 emission.

Visualizing excited state density differences with VMD

This task is almost identical to the visualization of the S0 → S1 excitation discussed above, but instead of

akdens ex

you need to call

akdens em

Blue color corresponds to volumes where the electron density increases during the S1/T1→ S0 emission.

Red color corresponds to volumes where the electron density decreases during the S1/T1→ S0 emission.

The plt from egrad file actually contains densities for a vertical excitation at the used geometry, but we can consider the vertical emission as the opposite case. Thus, the command "akdens em" simply reverses the sign of the blue and red isosurfaces in comparison to "akdens ex". This way, the interpretation of the colors remains the same (blue = increase, red = decrease).

  • No labels