How to calculate the Equilibrium Structure of a diatomic molecule

 In this tutorial, we will calculate the equilibrium structure of the Cl2 molecule.

The Cl2 molecule has only 2 atoms, and the structure is completely determined by the Cl−Cl bond length. Therefore we can determine the equilibrium structure by calculating the total energy as a function of the Cl−Cl distance.

The first step is to find a suitable pseudopotential for Cl. We go to http://www.quantum-espresso.org/pseudopotentials and look for Cl. We recognize LDA psedupotential by the label pz in the filename. Let us go for the following:

$ wget http://www.quantum-espresso.org/wp-content/uploads/upf_files/Cl.pz-bhs.UPF

 

We use the following input file of the Cl2 molecule to do the calculation:

cl2.in

&control
calculation = 'scf'
prefix = 'Cl2',
pseudo_dir = './',
outdir = './'
/
&system
ibrav = 1,
celldm(1) = 20.0,
nat = 2,
ntyp = 1,
ecutwfc = 100,
/
&electrons
conv_thr = 1.0d-8
/
ATOMIC_SPECIES
Cl 1.0 Cl.pz-bhs.UPF
ATOMIC_POSITIONS bohr
Cl 0.00 0.00 0.00
Cl 2.00 0.00 0.00
K_POINTS gamma

 Using ibrav = 1 we select a simulation box which is simple cubic, with lattice parameter celldm(1). Here we are choosing a cubic box of side 20 bohr (1 bohr = 0.529167 ˚A). The keyword gamma means that we will be sampling the Brillouin zone at the Γ point, that is k = 0. This is fine since we want to study a molecule, not an extended crystal. Note that we increased the planewaves cutoff to 100 Ry: this number was obtained from separate convergence tests.

We begin with a scf calculation to get the convergence.

$ mpirun -np 4 pw.x -i cl2.in |tee cl2.out  # for parallel calculation

$ pw.x -i cl2.in |tee cl2.out  # for serial calculation

 

Incidentally, from the output file of this run (say cl2.out) we can see the various steps of the DFT self-consistent cycle (SCF). For example if we look for the total energy: 

$ grep "total energy" cl2.out

     total energy              =     -55.58004522 Ry
     total energy              =     -55.80464971 Ry
     total energy              =     -55.84083719 Ry
     total energy              =     -55.84152314 Ry
     total energy              =     -55.84159990 Ry
     total energy              =     -55.84161846 Ry
     total energy              =     -55.84162026 Ry
!    total energy              =     -55.84162115 Ry
     The total energy is the sum of the following terms:

Here we see that the energy reaches its minimum in 8 iterations. The iterative procedure stops when the energy difference between two successive iterations is smaller than conv thr = 1.0d-8.

Now we calculate the total energy as a function of the Cl−Cl bond length. In the reference frame chosen for the input file above we have one Cl atom at (0,0,0) and one at (2,0,0) in units of bohr. Therefore we can vary the bond length by simply displacing the second atom along the x axis. In order to automate the procedure we can use the following script:

opt-bondlength.sh

#!/bin/sh

> cl2-bondlength.dat
echo -n "# " >> cl2-bondlength.dat
echo -n " Bondlength" >> cl2-bondlength.dat
echo  "  Energy" >> cl2-bondlength.dat


sed "s/2.00/NEW/g" cl2.in > tmp
for DIST in 2.2 2.4 2.6 2.8 3.0 3.2 3.4 3.6 3.8 4.0 4.2 4.4 4.6 
do
sed "s/NEW/${DIST}/g" tmp > cl2_${DIST}.in

mpirun -np 4 pw.x -i cl2_${DIST}.in |tee cl2_${DIST}.out  # for parallel calculation

# pw.x -i cl2_${DIST}.in |tee cl2_${DIST}.out  # for parallel calculation


echo -n ${DIST} >> cl2-bondlength.dat 
echo -n " " >> cl2-bondlength.dat

grep ! cl2_${DIST}.out | awk '{print $5}' >> cl2-bondlength.dat

done

We get the file cl2-bondlength.dat which contains the following data

#  Bondlength  Energy
2.2 -57.51390370
2.4 -58.54440037
2.6 -59.17256643
2.8 -59.55021751
3.0 -59.77201221
3.2 -59.89671133
3.4 -59.96077584
3.6 -59.98691108
3.8 -59.98945934
4.0 -59.97764683
4.2 -59.95749243
4.4 -59.93293867
4.6 -59.90658487

To the plot the data, we use the matplotlib package which is developed with python language. To install the package we run the following command:

sudo apt-get install python3-matplotlib

To the plot the file cl2-bondlength.dat we use the script Cl2-opt-bondlength.py

#!/usr/bin/env python3

import matplotlib.pyplot as plt
from matplotlib import rcParamsDefault
import numpy as np
#%matplotlib inline
plt.rcParams["figure.dpi"]=150
plt.rcParams["figure.facecolor"]="white"


file_name = input("Enter the file name:")
x, y = np.loadtxt(file_name, delimiter=' ', unpack=True)
plt.plot(x, y, "o-", markersize=5, label='Etot vs bond length')
plt.xlabel('bond length (Bohr)')
plt.ylabel('Etot (Ry)')
plt.legend(frameon=False)
plt.show()

To execute the script we run the following commands:

$ chmod +x Cl2-opt-bondlength.py

$ python3 Cl2-opt-bondlength.py

We get the following picture:

By zooming in the figure we find that the bond length at the minimum is 3.725 bohr = 1.97 ˚A. This value is 1.5% below the measured bond length of 1.99 ˚A.