How to Run Calculations for Simple Molecules

by Pasquale Pavone for exciting beryllium

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.

### 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 Tutorial scripts and environment variables. 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 = 0m56s

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.152852633
is =    2 (O), ia =    2 :    2.152852633

Distance between is =    2 (O), ia =    1 and
is =    1 (C), ia =    1 :    2.152852633
is =    2 (O), ia =    1 :    0.000000000
is =    2 (O), ia =    2 :    4.305705267

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

Therefore, the calculated bond lengths are

• dCO = 2.153 Bohr,
• dOO = 4.306 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 analize 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.