3. Tutorial for geometry optimisation

This tutorial will guide you step by step through a typical ONIOM calculation done with fromage. This page is meant to be self contained but for a refresher on the abbreviations used, visit the glossary. Everything discussed herein is elaborated upon in the rest of the documentation. The example calculation is rather computationally costly and time intensive. This is because we choose an model system with a distinctive excited state minimum and conical intersection geometry.

The example for this tutorial is the investigation of the emission properties of the (2-hydroxyphenyl)propenone (HP) crystal using the Ewald Embedded Cluster model (EEC) and Gaussian. The extension of this protocol to other programs is straightforward. The following steps will involve calculating the geometries and associated energies of the ground state, first excited state and S1 - S0 Minimal Energy Conical Intersection (MECI).

The necessary input files are in the tutorial directory. We follow the ONIOM nomenclature for the different regions of the embedded cluster. In other words the whole cluster is the “real system” and the central molecule is called the “model system”. We abbreviate the principal calculations involved in the ONIOM method as mh for the model system at a high level of theory, ml for the model system at a low level of theory and rl for the real system at a low level of theory. For the location of conical intersections, a fourth calculation is typically involved, mg, which indicates the model system at the high level of theory but in the ground state.

3.1. Calculation set up

The initial step is to generate the Ewald charge background as well as all of the necessary files needed for geometry optimisation.

3.1.1. Input

Make a directory where the preparatory calculation will take place:

mkdir opt_1/
cd opt_1/
cp path/to/fromage/tutorial/* .

The input files that you have just copied are:

  • cell.xyz

    A file containing the optimised positions of the atoms in the unit cell

  • config

    The fromage configuration file, described below

  • high_pop.log

    A Gaussian output file of a population analysis for one molecule using the high level of theory (in this case PBE 6-31G*)

  • low_pop.log

    The same but for the low level of theory (here, HF STO-3G)

  • *.template

    The template files for your upcoming ONIOM calculation. This is includes mh.template, ml.template, rl.template and mg.template. Note that the checkpoint file name is gck.chk in all cases and that the levels of theory match those of the population log files

3.1.2. Execution

If your installation was successful, all of the fromage scripts should be in your system path already. In that case, running the program simply involves typing:

fro_prep_calc.py

3.1.3. Output

After a few minutes, you will be greeted with a series of outputs:

  • prep.out

    Output file with some information about the setup calculation

  • mol.init.xyz

    The initial position of the model system

  • shell.xyz

    The molecules surrounding the model system

  • mh/ ml/ rl/ mg/

    Directories containing a .temp file each. For example mh/ contains mh.temp

  • ewald/

    The directory where the ewald calculation is run. The outputs in here are not important for this tutorial

3.2. Geometry optimisation

We will calculate the geometries and associated energies of the ground state minimum, first excited state minimum and S1 - S0 Minimal Energy Conical Intersection (MECI).

3.2.1. Ground state

3.2.1.1. Input

These are all the files needed for the geometry optimisation. Most of them were already generated from the previous step.

  • fromage.in (not required for this tutorial)

    The input file which contains the specifications for the geometry optimisation. This tutorial uses default specifications, so this is not required

  • mol.init.xyz

    See above

  • shell.xyz

    See above

  • mh/ ml/ rl/

    Directories containing their respective .temp files

3.2.1.2. Execution

An important part of calculations in fromage is the assignment of memory to each component calculation. Some times, depending on the system size and the combination of methods used, rl will need more memory than mh. Make sure to adapt the memory requested in all three .temp files to match the capacity of your system.

When this is ready, submit your job with the command:

fro_run.py

On the command line or in your job queue.

3.2.1.3. Output

You can expect this calculation to take a few hours depending on your computational resources. The convergence criterion of the optimisation is very strict by default so it is up to the user’s judgement whether they wish to abort the calculation once they have achieved a satisfactory precision.

  • fromage.out

    The main output file. This contains information about the energies and gradients at each step of the optimisation

  • geom_mol.xyz

    The positions of the model system throughout the optimisation

  • geom_clust.xyz

    The position of the real system throughout the optimisation. Only the model system will change

geom_mol.xyz should show very slight rearrangement of the molecule since its Gaussian-optimised ground state geometry is close its crystal.

3.2.1.4. Vertical excitation

To calculate the vertical excitation, make a new directory called exci/ and copy the mh.com file to it (it should now contain the optimised geometry):

mkdir exci/
cp mh/mh.com exci/
cd exci/

Then edit the mh.com file to remove the force keyword and add td(nstates=5,root=1). Now run Gaussian (or use a submission script):

g16 mh.com

This will give the excitation energies of the first five excited states, easily accessible with a judicious grep:

grep 'Excited State' mh.log

The output will look like this:

Excited State   1:      Singlet-?Sym    4.1338 eV  299.93 nm  f=0.6373  <S**2>=0.000
Excited State   2:      Singlet-?Sym    4.3687 eV  283.80 nm  f=0.0068  <S**2>=0.000
Excited State   3:      Singlet-?Sym    4.5466 eV  272.70 nm  f=0.0760  <S**2>=0.000
Excited State   4:      Singlet-?Sym    5.4041 eV  229.43 nm  f=0.0318  <S**2>=0.000
Excited State   5:      Singlet-?Sym    5.8322 eV  212.59 nm  f=0.0431  <S**2>=0.000

3.2.2. First excited state

We now wish to optimise the geometry of the molecule in the first excited state. The procedure will be almost identical to the one for the ground state optimisation but with an added keyword to mh.temp.

3.2.2.1. Input

First, copy the whole opt_1 directory to conserve the ground state data. Presuming you are still in opt_1/exci/, just type:

cd ../../
cp -r opt_1/ opt_2/
cd opt_2/

Now edit the file mh/mh.temp to add the keyword td(nstates=1,root=1).

And edit your mol.init.xyz to match the last geometry in geom_mol.xyz from your opt_1/ directory.

3.2.2.2. Execution

As usual, type:

fro_run.py

Or submit it to your job scheduler.

3.2.2.3. Output

This should typically take longer than your ground state calculation if you have succeeded in setting it up in a way that the limiting calculation is mh.

As described above, you will receive fromage.out, geom_mol.xyz and geom_clust.xyz.

This time, you should be able to see the excited state proton transfer in geom_mol.xyz as the optimised structure is in keto form.

3.2.3. MECI

3.2.3.1. Input

One final time, copy the whole directory:

cd ..
cp -r opt_2/ opt_3/
cd opt_3/

Edit the mol.init.xyz file with the final geometry of opt_2/geom_mol.xyz.

And in fromage.in, add a line at the bottom bool_ci 1. This turns on MECI search. Keep in mind that this calculation will use mg so change the memory requested in all of your .temp files accordingly.

3.2.3.2. Execution

Again:

fro_run.py

And wait however long it takes.

3.2.3.3. Output

The usual fromage.out, geom_mol.xyz and geom_clust.xyz will be generated.

fromage.out will contain different information, pertaining to the value and gradients of the penalty function which is being minimised instead of the energy.