Computational Chemistry, orca, programs

Calculations with Orca

This is a very quick post on how to run a calculation in Orca. In this case, it is an optimization of the iridium pentachloronitrosyl anion. This was a post that I wrote for my research group to encourage the use of Orca for quantum mechanical calculations. I will write this up further in due course.

For more detailed information, check out the Orca manual and the Orca Input Library.

Continue reading

Computational Chemistry

Driver calculations in Orca

In this example, we will set up and run a driver calculation in Orca. In this type of calculation, one variable (an angle or bond length) is varied while the energy is minimised. This type of calculation can be useful for exploring reaction pathways. In this example, we will use the dissociation of a carbonyl ligand from a rhodium-phosphine complex. These types of complexes are common in industrial chemistry for a variety of catalytic reactions; however, typically more complicated phosphine or phosphite ligands, either mono- or bidentate, are used. For computational simplicity, I am using phosphine (PH3) as the ligand.

Setting up the calculation

For simplicity, we will use a split valence basis on all atoms except rhodium. On rhodium, we will use a Stuttgart-Dresden ECP with the def2-tzvp basis set modelling the valence electrons. The calculation uses BP86 and Grimme’s D3 dispersion correction with Becke-Johnson damping. Below is the input for this driver calculation.

! noautostart slowconv tightscf d3bj opt bp86 pal2 def2-svp def2-svp/j ri ecp{def2-sd=[Rh], def2-tzvp, def2-tzvp/j}
%geom scan
B 0 11 = 1.909, 4.0, 8
*xyz 0 1

Concerning the topic of this post, i.e., the driver calculation, the most important part of this input file is the part describing the bond length that we are ‘driving’. In this example, the Rh-C bond length. In the XYZ coordinates, Rh and C (of carbonyl) are the first and twelfth atoms, respectively. Orca numbers atoms from 0 (zero) upwards, so it is important to remember that the nth numbered atom is the n-1th atom in Orca. Here, the command ‘B 0 11 = 1.909, 4.0, 8’ indicates that a bond (B) is changed from 1.909 to 4 Å in eight steps, i.e., 2.6 Å per step.

In addition, the ‘simple’ command line contains several additional elements. The first command ‘noautostart’ indicates that we do not want to read from an existing wavefunction file, should it exist. ‘slowconv’ and ‘tightscf’ are commands that control the SCF calculation, and these are useful for transition metal complexes. The command ‘d3bj’ indicates that we want to model dispersion interactions using Grimme’s D3 correction with Becke-Johnson damping, and ‘opt’ indicates that this is an optimisation, that is, optimising everything other than the coordinate being driven. ‘pal2’ indicates that this is a parallel calculation, running on two openmpi processes. ‘def2-svp def2-svp/j ri ecp{def2-sd=[Rh], def2-tzvp, def2-tzvp/j}’ are commands related to the basis sets on the atoms in this system. A split basis (def2-svp) is being used on all atoms; however, on Rhodium, an electron core potential is being used to model the core electrons, and a triple-zeta basis is used for the valence electrons. In addition, the RI approximation is used to speed up the calculations, and this si the reason for the basis sets appended with ‘/j’, i.e., these are the density fitting basis sets.

About the calculation

The GIF below shows the eight steps of the relaxed driver calculation; for fun, it is a stereogram. The Rh-C bond was lengthened in steps of about 0.3 A. At each step, an optimization was carried out, but the Rh-C distance was kept constant.


As the carbonyl ligand moves away from the rhodium atom, the Rh-C-O bond angle changes. Starting off linear at 177.2 degrees; it becomes bent with an angle of about 109.

Charting the Rh-C-O bond angle versus the Rh-C bond length.

Of course, after a certain point there is no interaction between Rh and the C of the carbonyl ligand; this is given by the dissociation limit of the complex. The Rh-C-O tilt angle shows how the bonding interaction changes from a ligand sigma bonded with π-backbonding to a free ligand, and this is nicely shown by the shortening of the C-O bond length from 1.16 A to 1.14 A as the carbonyl ligand dissociates.

Analysing the data

Energy profile for the driver calculation.

The results of the calculation are shown to the left. As you can see, the CO-metal bond elongates, followed by a tilting of CO and its slow drift away from the complex. The energy of the complex increases until it reaches the dissociation limit; here, the energy levels off to a constant value.

In fact, we can fit the energy profile to the Morse potential and get an idea of the dissociation limit and the equilibrium distance for the Rh–C bond length.

The Morse equation is y = D \times(1-exp(-a\times (x-e)))^2 where D is the dissociation limit, a is the well depth, and e is the equilibrium distance.

Fitting of a Morse curve to the driver data.
Fitting of a Morse curve to the driver data.

Fitting this equation to the data, The morse curve gives values for D, a, and e of 26.92 kcal/mol, 2.81 kcal/mol, and 1.92 A, respectively.

chemistry, Computational Chemistry

Checking on geometry optimisations in Gaussian.

This is a useful oneliner that uses egrep to check on the status of a geometry optimisation. It is taken from Exploring chemistry with electronic structure methods: a guide to using Gaussian*, a very useful, albeit old, book about Gaussian 98. I usually copy the code below into a bash script, change its mode with chmod** and then place it in my personal script folder (which is indicated in my path). You could alternatively just put the egrep line as an alias in your shell profile.

egrep 'out of| SCF Don|Converged| NO |YES|exceeded' $file | grep -v '\\\\'

*Foresman, James B., and Æleen Frisch. Exploring chemistry with electronic structure methods: a guide to using Gaussian. Gaussian, 1993.

** chmod u+x

chemistry, Computational Chemistry, pymol

Pymol for computational chemistry.

A picture is worth a thousand words.

So goes the saying, and it is true.

The pymol wiki has lots of ideas for styling atoms. Below are instructions for styling atoms, as in the image shown below, which I think produces an image that is particularly clear (and pretty) and suited to publication in journal articles.

hide everything
show sticks
show spheres
set stick_radius, .07
set sphere_scale, .18
set sphere_scale, .13, elem H
set bg_rgb=[1, 1, 1]
set stick_quality, 50
set sphere_quality, 4
color gray85, elem C
color red, elem O
color slate, elem N
color gray98, elem H
set stick_color, black
set ray_trace_mode, 1
set ray_texture, 2
set antialias, 3
set ambient, 0.5
set spec_count, 5
set shininess, 50
set specular, 1
set reflect, .1
set dash_gap, 0
set dash_color, black
set dash_gap, .15
set dash_length, .05
set dash_round_ends, 0
set dash_radius, .05

This can be copied directly into the pymol command line, followed by enter. If you have non-bonded atoms such as a central metal with attached ligands, for which pymol sometimes does not detect, you can use the measurement tool and then remove the labels and change dash_radius to a very small number; alternatively, use the bond command.

An image of xantphos rendered in PyMol
An image of xantphos rendered in PyMol and in stereo. Cross your eyes to see a three-dimensional image.

Another nice effect is to change the display types of certain atoms. For example, I have highlighted the  coordination of the oxygen atoms in the crown ether molecules by hiding the spheres:

 hide spheres, elem C

Also, shown below is a GIF animation an animation of an optimization of the crown ether complex.