How to Run Calculations for Simple Molecules

by Pasquale Pavone & Caterina Cocchi for exciting boron

Purpose: In this tutorial you will learn how to set up and execute exciting calculations for simple molecules such as the carbon dioxide (CO2), water (H2O), and pyridine (C5H5N) molecule. It is also shown how to calculate and visualize molecular orbitals, with explicit application to the HOMO and LUMO of carbon dioxide.

### 0. Define relevant environment variables

Read the following paragraphs before starting with the rest of this tutorial!

Before starting, be sure that relevant environment variables are already defined as specified in How to set environment variables for tutorials scripts. Here is a list of the scripts which are relevant for this tutorial with a short description.

• EXECUTE-single.sh: (Bash) shell script for running a single exciting calculation.
• PLOT-relaxdistance.py: Python visualization tool for following relative atomic coordinates of atoms during the relaxation process.
• PLOT-status.py: Python visualization tool for RMS deviations of the SCF potential as a function of the iteration number during the SCF loop.
• PLOT-maxforce.py: Python visualization tool for the maximum amplitude of the force on the atoms during relaxation.
• PLOT-torque.py: Python visualization tool for the total torque during relaxation.
• PLOT-centerofmass.py: Python visualization tool for the cartesian components of the position of the center of mass during relaxation.

From now on the symbol \$ will indicate the shell prompt.

Requirements: Bash shell. Python numpy, lxml, matplotlib.pyplot, and sys libraries.

### 1. Carbon-dioxide molecule

The CO2 molecule is a linear one with the two oxygen atoms located symmetrically with respect to the central carbon atom, as can be seen in the following figure taken from the wikipedia page of Carbon dioxide. The experimental bond lengths are

• dCO = 2.198 Bohr,
• dOO = 4.396 Bohr.
##### i) Isolated molecules are treated as special periodic systems

In order to perform the calculation, we have to let exciting know that we would like to consider an isolated molecule with a certain atomic configuration. However, all calculations performed by exciting assume a three-dimensional periodicity along the directions defined by the lattice vectors. Therefore, when considering a molecule one has to prepare a special three-dimensional crystal where the distances between a molecule and each of its replicas are large enough so that each molecule behaves as it would be isolated.

The procedure one has to follow in this case is illustrated in the next picture in the simplification of only two dimensions. The first thing to do is to determine the size (expressed by dx, dy, dz along the three cartessian directions) of the molecule you are investigating, as illustrated in the left panel of the picture above. Next, some vacuum (expressed by vx, vy, vz) is added, as you can see above in the middle panel. The amount of vacuum has to be chosen in such a way to prevent intermolecular interactions. Finally, a cartesian periodic cell is created with basis vectors of amplitude (dx+vx, dy+vy, dz+vz) along the corresponding cartesian axis (see right panel).

##### ii) Preparation of the input file

The first step is to create a directory for each system that you want to investigate. Thus, we will create a directory CO2 and we move inside it.

``````\$ mkdir CO2
\$ cd CO2```
```

An example for an exciting input file (input.xml) which is needed for the calculation of the optimized structure of the CO2 molecule is given in the following.

```<input>

<title>CO2 molecule</title>

<structure speciespath="\$EXCITINGROOT/species" cartesian="true">
<crystal>
<basevect> 10.0    0.0    0.0 </basevect>
<basevect>  0.0    5.0    0.0 </basevect>
<basevect>  0.0    0.0    5.0 </basevect>
</crystal>
<species speciesfile="C.xml" rmt="1.0">
<atom coord=" 0.00  0.00  0.00" lockxyz="true true true"/>
</species>
<species speciesfile="O.xml" rmt="1.0">
<atom coord=" 2.50  0.00  0.00" lockxyz="false true true"/>
<atom coord="-2.50  0.00  0.00" lockxyz="false true true"/>
</species>
</structure>

<groundstate
xctype="LDA_PW"
gmaxvr="16"
rgkmax="4"
outputlevel="low">
</groundstate>

<relax
printtorque="true"
history="true">
</relax>

</input>
```

Some of the elements and attributes appearing in this input file should be already familiar to you from other tutorials. Here, we will discuss in details only those elements and attributes which are significant in molecule calculations.

First, we analyze the first line of the element structure.

```   <structure ... cartesian="true">
```

Here, the new attribute cartesian appears. If cartesian is set to "true" the atomic coordinates which are specified by lines like

```         <atom coord=" 2.50  0.00  0.00" .../>
```

are assumed to be given in cartesian coordinates. This option is particularly useful when dealing with molecules.

Furthermore, the basis vectors, which define the periodic cell including the molecule and the vacuum, are taken in the most simple way to lay along the cartesian axis. The amplitude of each basis vector is chosen in such a way to accommodate the size of the molecule in its direction and the amount of vacuum which we desire to add. The basis vectors shown in our example are

```         <basevect> 10.0    0.0    0.0 </basevect>
<basevect>  0.0    5.0    0.0 </basevect>
<basevect>  0.0    0.0    5.0 </basevect>
```

and correspond to a linear molecule of size of about 5 Bohr lying along the x axis and a vacuum of 5 Bohr along each cartesian direction.

The next new feature that can be seen in the input file is the attribute lockxyz which is always located inside the atom element.

```         <atom coord=" 2.50  0.00  0.00" lockxyz="false true true"/>
```

The optional attribute lockxyz is connected to the optimization of the internal degree of freedom of a system (molecule or crystal) which is performed when using the element relax. In particular, if for one atom one of the three logical values which are given to lockxyz is set to true, this means that the corresponding cartesian component of the considered atom is kept fixed during the optimization process. As an example, writing

```         <atom coord=" 2.50  2.50  2.50" lockxyz="false false true"/>
```

means that, for this atom, the x and y cartesian coordinates are allowed to change during the relaxation, while the z component of the position vector is kept fixed. Furthermore, when considereing molecules, it is always convenient to keep fixed the position of one of the atoms of the molecule in order to avoid the translational motion of the molecule as a whole. We impose this condition on the central carbon atom.

```      <species speciesfile="C.xml" ...>
<atom coord=" 0.00  0.00  0.00" lockxyz="true true true"/>
</species>
```

The next important setting, visible in the input file input.xml is that the muffin-tin radii for the atomic species are given explicitly and assume values which are significantly smaller than the default values.

```...
<species speciesfile="C.xml" rmt="1.0">
...
<species speciesfile="O.xml" rmt="1.0">
...
```

This is necessary due to the short bond lengths between the atoms in the molecule.

Notice that, contrary to all the crystal structures you have considered in these tutorials, the groundstate element does not contain explicitly the attribute ngridk. This is completely equivalent to set ngridk to "1 1 1" as in the following example.

```   <groundstate
ngridk="1 1 1"
...
</groundstate>
```

This choice is necessary for dealing with a molecule for which one assumes that no interaction is present between the molecule and its periodic replicas.

The next thing that can be noticed in the input file, is that the attribute rgkmax of the groundstate element has the value "4".

```   <groundstate
...
rgkmax="4"
...
</groundstate>
```

This value is significantly smaller than typical values used for crystals (where converged calculations require a value of rgkmax around "7" or "8"). Due to presence of the vacuum region which reduce the number of plane waves necessary to describe the interstitial region in the LAPW method, setting rgkmax to "4" (or even to smaller values) is typical for molecule calculations.

The final remark concerns the attribute printtorque (of the element relax) which has been set to "true".

```   <relax
printtorque="true"
...
</relax>
```

At each optimization step, this setting produces in the output file "INFO.OUT" two lines which contain the cartesian components of the center of mass and total torque, respectively. These lines should look like the following.

``````...
Center of mass  :     0.00001466   -0.00050662    0.00000000  (cartesian)
Total torque    :     0.00000000    0.00000000    0.00000798  (cartesian)
...```
```

The center of mass and the total torque are well defined only when dealing with isolated systems such as those considered in this Tutorial. In particular, towards the end of the optimization run, non vanishing values for the components of the total torque (numerically, larger than about 10-4 Ha) may indicate that the molecule as a whole has started to rotate or oscillate. This would prevent the relaxation to reduce the atomic forces below the desidered threshold.
##### iii) Execute calculations

N.B.: Do not forget to replace in the input.xml the string "\$EXCITINGROOT" by the actual value of the environment variable \$EXCITINGROOT using the command

````\$ SETUP-excitingroot.sh`
```

To perform the struture optimization of the CO2 molecule you have to run the script EXECUTE-single.sh in the directory where the input file input.xml has been created. If you want to save the results of this run in the subdirectory 1-CO2 you should run as in the following.

``````\$ EXECUTE-single.sh 1-CO2

===> Output directory is "1-CO2" <===

Running exciting for file input.xml -------------------------------------

Elapsed time = 1m12s

Run completed for file input.xml ----------------------------------------

\$```
```
##### iv) Analyze results

The analysis of the result of the optimization run can be visualized with the help of the scripts already introduced in Simple examples of structure optimization.

The maximum value of the force acting on the atoms can be visualized using the script PLOT-maxforce.py

````\$ PLOT-maxforce.py 1-CO2`
```

The output of this script (files PLOT.ps or PLOT.png) should look like the following. We can see that the optimization required 5 steps to reach the target minimum force.

Further information is provided using the script PLOT-relaxdistance.py.

````\$ PLOT-relaxdistance.py 1-CO2`
```

In this case, the output reproduces the distance between one of the oxygen atoms and the central carbon and looks like Notice that in this case crystal axes are parallel to the cartesian directions. One can also notice that, as expected, the y (2) and z (3) components do not change during the optimization, while the bond length along the x (1) axis is reduced from the initial value and reaches a value slightly larger than 2.1 Bohr.

The complete information about the optimized bond lengths is found in the file BONDLENGTH_OPT.OUT in the 1-CO2 directory. Typing

````\$ cat 1-CO2/BONDLENGTH_OPT.OUT`
```

you should have an output on the screen similar to

``````Distance between is =    1 (C), ia =    1 and
is =    1 (C), ia =    1 :    0.000000000
is =    2 (O), ia =    1 :    2.152739010
is =    2 (O), ia =    2 :    2.152739010

Distance between is =    2 (O), ia =    1 and
is =    1 (C), ia =    1 :    2.152739010
is =    2 (O), ia =    1 :    0.000000000
is =    2 (O), ia =    2 :    4.305478021

Distance between is =    2 (O), ia =    2 and
is =    1 (C), ia =    1 :    2.152739010
is =    2 (O), ia =    1 :    4.305478021
is =    2 (O), ia =    2 :    0.000000000```
```

Therefore, the calculated bond lengths are

• dCO = 2.153 Bohr,
• dOO = 4.305 Bohr.
##### Exercises
• Perform the optimization of the CO2 molecule starting from a planar configuration by setting initial coordinates and the attribute lockxyz of the two Oxygen atoms to the following (Notice: in comparison to the previous input.xml, the lock on the y component has been deleted!).
```      <species speciesfile="O.xml" rmt="1.0">
<atom coord=" 2.50  0.50  0.00" lockxyz="false false true"/>
<atom coord="-2.50  0.50  0.00" lockxyz="false false true"/>
</species>
```
• Investigate the dependence of the relaxed structure on the amount of vacuum used in the calculation. Try both increasing and decreasing the value used in the example.
• Investigate the dependence of the relaxed structure on the value of the attribute rgkmax.
• Change the initial configuration (using low symmetry displacements) as in the following and perform again the optimization in a directory named CO2-low-symmetry.
```      <species speciesfile="O.xml" rmt="1.0">
<atom coord=" 2.45  0.60  0.00" lockxyz="false false true"/>
<atom coord="-2.50  0.40  0.00" lockxyz="false false true"/>
</species>
```
• Concerning the results in CO2-low-symmetry. How long does it take the optimization run in this case? Why?
• Have a look to the behaviour of the total torque during the optimization run using the script PLOT-torque.py. If results of this optimization run are in the directory CO2-low-symmetry write
````\$ PLOT-torque.py CO2-low-symmetry`
```
• Have a look to the behaviour of the position of the center of mass during the optimization run using the script PLOT-centerofmass.py. If results of this optimization run are in the directory CO2-low-symmetry write
````\$ PLOT-centerofmass.py CO2-low-symmetry`
```

### 2. Water molecule

The H2O molecule is a planar one with the two hydrogen atoms located symmetrically with respect to the central oxygen atom, as can be seen in the following figure taken from the wikipedia page of water. The experimental bond lengths are

• dHO = 1.811 Bohr,
• dHH = 2.863 Bohr.

The experimental value of the HOH angle is

• = 104.45 degrees

Also in this case , the first step is to create a new directory for the system that you want to investigate. Thus, we will create a directory H2O and we move inside it.

``````\$ mkdir H2O
\$ cd H2O```
```

An example of input file (input.xml) for the optimization of the atomic structure of the water molecule is the following.

```<input>

<title>Water molecule</title>

<structure speciespath="\$EXCITINGROOT/species" cartesian="true">
<crystal>
<basevect> 10.0    0.0    0.0 </basevect>
<basevect>  0.0   10.0    0.0 </basevect>
<basevect>  0.0    0.0    5.0 </basevect>
</crystal>
<species speciesfile="O.xml" rmt="1.0">
<atom coord="0.00  0.00  0.00" lockxyz="true  true  true"/>
</species>
<species speciesfile="H.xml" rmt="0.7">
<atom coord="2.50  0.50  0.00" lockxyz="false false true"/>
<atom coord="0.50  2.50  0.00" lockxyz="false false true"/>
</species>
</structure>

<groundstate
xctype="LDA_PW"
gmaxvr="16"
rgkmax="3"
outputlevel="low">
</groundstate>

<relax
taubfgs="0.5"
printtorque="true"
history="true">
</relax>

</input>
```

Replace in the input.xml the string "\$EXCITINGROOT" by the actual value of the corresponding environment variable.

##### Exercises
• Try to analyze this input file in the light of the experience gained with the previous molecule.
• Perfom the relaxation calculation.
• Increase and decrease the vacuum around the molecule (try each direction independently), what happens? (Hint: use the visualization tool PLOT-status.py in order to follow the convergence towards the targets for the SCF cycles.)
• Change the initial configuration (allowing for low symmetry displacements) and perform again the optimization.
• Use different values of the attribute method of the element relax for the structure optimization.
• Perform the convergence study of the optimized structure as a function of the value of the parameter rgkmax.

### 3. Pyridine molecule

Pyridine is a basic organic compound with the chemical formula C5H5N with an atomic structure similar to benzene, as can be seen in the following figure taken from the wikipedia page of Pyridine. Create a directory Pyridine and move inside it. An example of input file (input.xml) for the optimization of the atomic structure of the pyridine molecule is the following.

```<input>

<title>Pyridine</title>

<structure speciespath="\$EXCITINGROOT/species" cartesian="true">
<crystal>
<basevect> 14.0      0.0     0.0 </basevect>
<basevect>  0.0     14.0     0.0 </basevect>
<basevect>  0.0      0.0     5.0 </basevect>
</crystal>
<species speciesfile="N.xml" rmt="1.05">
<atom coord=" 0.00    0.00    0.00" />
</species>
<species speciesfile="C.xml" rmt="1.05">
<atom coord="-2.18    1.30    0.00" />
<atom coord="-2.25    3.99    0.00" />
<atom coord=" 0.00    5.25    0.00" />
<atom coord=" 2.25    3.99    0.00" />
<atom coord=" 2.18    1.30    0.00" />
</species>
<species speciesfile="H.xml" rmt="0.9">
<atom coord="-3.93    0.10    0.00" />
<atom coord="-4.03    5.27    0.00" />
<atom coord=" 0.00    7.67    0.00" />
<atom coord=" 4.03    5.27    0.00" />
<atom coord=" 3.93    0.10    0.00" />
</species>
</structure>

<groundstate
xctype="LDA_PW"
gmaxvr="16"
rgkmax="3.0">
</groundstate>

<relax
printtorque="true"
history="true">
</relax>

</input>
```

Replace in the input.xml the string "\$EXCITINGROOT" by the actual value of the corresponding environment variable.

##### Exercises
• Perform similar investigations as for the other molecules.

### 4. Visualization of molecular orbitals of CO2

In this section, we will learn how to plot molecular orbitals (MOs). More specific details about visualization tools available in exciting are reported in How to visualize Kohn-Sham states. As a test example, we will consider here the CO2 molecule. Hence, we can refer to the directory created for the runs described in Section 1.

In order to create a 3D plot of the MO isosurfaces, we need to perform a ground-state calculation, which includes the elements properties and wfplot (see How to visualize Kohn-Sham states). The input file we are going to use looks like the following:

```<input>

<title>CO2 molecule: Plot example</title>

<structure speciespath="\$EXCITINGROOT/species" cartesian="true">
<crystal>
<basevect> 10.0    0.0    0.0 </basevect>
<basevect>  0.0    5.0    0.0 </basevect>
<basevect>  0.0    0.0    5.0 </basevect>
</crystal>
<species speciesfile="C.xml" rmt="1.0000">
<atom coord=" 0.0000000000   0.0000000000   0.0000000000" />
</species>
<species speciesfile="O.xml" rmt="1.0000">
<atom coord=" 2.1527390104   0.0000000000   0.0000000000" />
<atom coord="-2.1527390104   0.0000000000   0.0000000000" />
</species>
</structure>

<groundstate
xctype="LDA_PW"
gmaxvr="16"
rgkmax="4"
outputlevel="low">
</groundstate>

<properties>
<wfplot>
<kstlist>
<pointstatepair>1 8</pointstatepair>
</kstlist>
<plot3d>
<box grid="40 40 40">
<origin coord="-0.5 -0.5 -0.5"/>
<point  coord=" 0.5 -0.5 -0.5"/>
<point  coord="-0.5  0.5 -0.5"/>
<point  coord="-0.5 -0.5  0.5"/>
</box>
</plot3d>
</wfplot>
</properties>

</input>
```

Please notice that in the structure element we replaced the atomic positions given in the input file in Section 1, with the optimized positions obtained from the relaxation performed in the same section. As usual, before running any calculation, replace in the input.xml the string "\$EXCITINGROOT" by the actual value of the corresponding environment variable.

In the input file shown above, the element wfplot inside properties triggers the calculation of the MO. wfplot includes the plot3D subelement, which in turn enables the 3D visualization of the orbital. The state which one wants to plot is defined by the two elements kstlist and pointstatepair. Inside the element pointstatepair a pair with k-point (first integer index) and eigenvalue (second integer index) number must be specified. As we are dealing here with an isolated molecule, only 1 k-point is present, therefore, the first index is always 1 for molecules. The second index, related to the eigenvalue labelling, is chosen as the highest occupied MO (HOMO), which corresponds to 8 in the CO2 molecule, as indicated in the file EIGVAL.OUT, where a map of eigenvalues is stored. Occupied (unoccupied) energy levels are identified by a nonvanishing (vanishing) number in the neighboring column. Kohn-Sham states with the same energy are degenerate. In particular, we notice that the HOMO of the CO2 molecule is doubly degenerate.

Additional elements are box, origin, and point. The element box is specified by the attribute grid, which defines the quality of the plot: The larger number of points in this mesh, the more resolved the resulting plot. Be aware that the increase of the number of grid points impacts on the computational costs! The elements origin and point allow to define, respectively, the origin and the vectors spanning the space region where the squared norm wave function is calculated. The values of the attributes coord are expressed in lattice coordinates and are chosen here in order to correctly place the orbital with respect to the atomic positions.

##### i) Visualization of the first HOMO

We now run the calculation as follows:

````\$ EXECUTE-single.sh HOMO_1`
```

When the run is completed, move to the newly created directory HOMO_1.

````\$ cd HOMO_1`
```

Inside this directory, we find the file WF3D.xml where the MO is stored. We are going to visualize it with XCrySDen. As a first step, we transform WF3D.xml into a xsf file, by running the following template:

````\$ xsltproc \$EXCITINGVISUAL/plot3d2xsf.xsl WF3D.xml > WF3D.xsf`
```

In order to visualize also the underlying structure of the molecule, we transform the atomic coordinates of the input file into a xsf file:

````\$ xsltproc \$EXCITINGCONVERT/xmlinput2xsf.xsl input.xml > input.xsf`
```

Finally, we merge the two files as follows:

````\$ cat input.xsf WF3D.xsf > HOMO_1.xsf`
```

Applying the procedure reported in How to visualize Kohn-Sham states for the visualization of xsf files using XCrySDen, we can now plot the isosurface of the probability density of the HOMO.

````\$ xcrysden --xsf HOMO_1.xsf &`
```

Using an isovalue one order of magnitude smaller than the maximum value, the following result is obtained: From this plot we observe that the HOMO of CO2 is correctly localized on the oxygen atoms of the molecule.

##### ii) Visualization of the second HOMO

As we noticed already from the file EIGVAL.OUT that the molecule CO2 has a doubly degenerate HOMO, we now plot the Kohn-Sham state corresponding to index 7 in the file EIGVAL.OUT. To do so, we first move to the parent directory.

````\$ cd ..`
```

Here, we modify inside the existing file input.xml the second index in the element pointstatepair as follows:

```...
<kstlist>
<pointstatepair>1 7</pointstatepair>
</kstlist>
...
```

Now we run again exciting using the command

````\$ EXECUTE-single.sh HOMO_2`
```

and, after the run is completed, move to the directory HOMO_2.

````\$ cd HOMO_2`
```

As a result of the calculation run, another file WF3D.xml is generated. Using the same procedure described above we can generate the file HOMO_2.xsf, which can be visualized using XCrySDen obtaining the following plot. The degeneracy of the two HOMO orbitals is evident: The shape of the isosurface is the same, while its orientation with respect to the coordinate system is rotated by 90 degrees.

##### iii) Visualization of the LUMO

Finally, we also calculate and plot the lowest unoccupied MO (LUMO) of CO2. The procedure is now known.

````\$ cd ..`
```

The LUMO is non degenerate and is labeled with the index 9 in the file EIGVAL.OUT. We modify the input file accordingly

```...
<kstlist>
<pointstatepair>1 9</pointstatepair>
</kstlist>
...
```

and run again exciting:

``````\$ EXECUTE-single.sh LUMO
\$ cd LUMO```
```

Repeating the standard visualization procedure inside the directory LUMO, we obtain the file LUMO.xsf which originates the following image (using as isovalue 0.015). ##### Exercises
• Inspecting the file EIGVAL.OUT, plot the probability density of deeper occupied and higher unoccupied MOs than HOMO and LUMO, respectively. How do high unoccupied MOs look like?
• Pay attention to the pairs of degenerate states: what is the character and the spatial distribution of these orbitals?
• Increase the amount of vacuum in the unit cell, by increasing the lattice vectors \$a\$, \$b\$, and \$c\$ to 12, 7, and 7 Bohr, respectively. Is there an influence on the MOs? What happens in particular to the empty states?
• What happens to the plots by choosing a different isovalue?
• Visualize the HOMO and LUMO for the water and pyridine molecule, too.