How to Run Calculations for Simple Molecules¶
by Pasquale Pavone & Caterina Cocchi for exciting sodium
(Jupyter notebook by Tejus Rohatgi, Mara Voiculescu & Martin Kuban)
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. Before Starting¶
Read the following paragraphs before starting with the rest of this tutorial!
Before running any Jupyter tutorials, please refer to the 00_before_starting.md document on how to correctly set up the environment. This only needs to be done once. After which, the venv can be (re)activated from exciting's root directory:
source $EXCITINGROOT/tools/excitingjupyter/venv/excitingvenv/bin/activate
Here is a list of the scripts which are relevant for this tutorial with a short description.
execute.single.py: Python script for running a single exciting calculation.convert_xml2xsf.py: Python script for converting xml to xsf files.sduring the relaxation process.plot.bondlength.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.
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 cartesian directions) of the molecule you are investigating, as illustrated in the left panel of the picture above. Next, some vacuum (expressed by $v_x$, $v_y$, $v_z$) 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 ($d_x+v_x$, $d_y+v_y$, $d_z+v_z$) 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.
%%bash
mkdir -p 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 detail 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 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 considering 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 in general 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-3 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 desired threshold.
%%bash
cd CO2
python3 -m excitingscripts.setup.excitingroot
cd ..
To perform the structure optimization of the CO2 molecule you have to run the script excitingscripts.setup.excitingroot.py 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 CO2-test-1 you should run as in the following.
%%bash
mkdir -p CO2/CO2-test-1
cd CO2/CO2-test-1
cp ../input.xml ./
python3 -m excitingscripts.execute.single -f input.xml
cd ../../
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 excitingscripts.plot.maxforce.py
%%bash
python3 -m excitingscripts.plot.maxforce -r CO2/CO2-test-1/
The output of this script (file PLOT.png) should look like the following.
We can see that the optimization required 4 steps to reach the target minimum force.
Further information is provided using the script excitingscripts.plot.bondlength.py.
%%bash
python3 -m excitingscripts.plot.bondlength CO2/CO2-test-1 1 2
In this case, the output reproduces the bond length between atom "1" (the central carbon atom) and atom "2" (the first oxygen atom as listed in input.xml) and looks like
Notice that in this case crystal axes are parallel to the cartesian directions. One can also notice that the bond length is reduced from the initial value and reaches a value slightly larger than 2.15 Bohr.
The complete information about the optimized bond lengths is found in the file BONDLENGTH_OPT.OUT in the CO2-test-1 directory. Typing
%%bash
cat CO2/CO2-test-1/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.152770757
is = 2 (O), ia = 2 : 2.152770757
Distance between is = 2 (O), ia = 1 and
is = 1 (C), ia = 1 : 2.152770757
is = 2 (O), ia = 1 : 0.000000000
is = 2 (O), ia = 2 : 4.305541513
Distance between is = 2 (O), ia = 2 and
is = 1 (C), ia = 1 : 2.152770757
is = 2 (O), ia = 1 : 4.305541513
is = 2 (O), ia = 2 : 0.000000000
Therefore, the calculated bond lengths are
- dCO = 2.153 Bohr,
- dOO = 4.306 Bohr.
- Perform the optimization of the CO2 molecule starting from a planar configuration by setting initial coordinates and the attribute
lockxyzof 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 the directory CO2-test-2.
<species speciesfile="O.xml" rmt="1.0">
<atom coord=" 2.45 0.05 0.00" lockxyz="false false true"/>
<atom coord="-2.55 -0.05 0.00" lockxyz="false false true"/>
</species>
Concerning the results in CO2-test-2. How long does it take to run the optimization in this case? Why?
Have a look to the behavior of the total torque during the optimization run using the script
excitingscripts.plot.torque.py. If results of this optimization run are in the directory CO2-test-2 writepython3 -m excitingscripts.plot.torque CO2-test-2By considering the total torque, compare the results of the runs CO2-test-1 and CO2-test-2. Why are they different?
Have a look to the behaviour of the position of the center of mass during the optimization run using the script
excitingscripts.plot.centerofmass.py. If results of this optimization run are in the directory CO2-test-2 writepython3 -m excitingscripts.plot.centerofmass CO2-test-2
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
$\theta$HOH= 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.
%%bash
mkdir -p 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
printtorque="true"
history="true">
</relax>
</input>
Replace in the input.xml the string $EXCITINGROOT by the actual value of the corresponding environment variable.
%%bash
cd H2O
python3 -m excitingscripts.setup.excitingroot
cd ..
Try to analyze this input file in light of the experience gained with the previous molecule.
Perform the relaxation calculation.
Increase and decrease the vacuum around the molecule (try each direction independently), what happens? (Hint: use the visualization tool
excitingscripts.plot.status.pyin 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
methodof the elementrelaxfor 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.
%%bash
mkdir -p Pyridine
Create a directory Pyridine. 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.
%%bash
cd Pyridine
python3 -m excitingscripts.setup.excitingroot
cd ..
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 the tutorial 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.
%%bash
mkdir -p CO2/HOMO_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.0">
<atom coord=" 0.00000 0.0 0.0" />
</species>
<species speciesfile="O.xml" rmt="1.0">
<atom coord=" 2.15277 0.0 0.0" />
<atom coord="-2.15277 0.0 0.0" />
</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 labeling, 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 in the number of grid points impacts the computational costs! The elements origin and point allow us 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.
%%bash
cd CO2/HOMO_1
python3 -m excitingscripts.setup.excitingroot
python3 -m excitingscripts.execute.single -f input.xml
cd ../../
When the run is completed, move to the newly created directory 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:
%%bash
cd CO2/HOMO_1
python3 -m excitingscripts.convert_xml2xsf -d 3D -f WF3D.xml
cd ../../
In order to visualize also the underlying structure of the molecule, we transform the atomic coordinates of the input file into a xsf file:
%%bash
cd CO2/HOMO_1
python3 -m excitingscripts.convert_xml2xsf -d 3D -f input.xml
cd ../../
Finally, we merge the two files as follows:
%%bash
cd CO2/HOMO_1
cat input.xsf WF3D.xsf > HOMO_1.xsf
cd ../../
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.
%%bash
cd CO2/HOMO_1
xcrysden --xsf HOMO_1.xsf >/dev/null 2>&1 &
cd ../../
$\Rightarrow$ See here the complete procedure to visualize 3D-plots in xsf format using XCrySDen
After running the previous command a XCrySDen window will appear on the screen. If you want to visualize the isosurface, use the following procedure:
- Press on the main window Tools -> Data Grid. A new window appears.
- Press OK: The previous window will disappear and a new window will appear with the control panel of the data grid to be plotted.
- In the middle of the left column a box appears, indicating Minimum grid value, Maximum grid value, and Isovalue next to an empty box. In this box the user has to type the isovalue with which the isosurface is to be plotted. As a rule of thumb, 10% of the Maximum grid value is a reasonable value. Type the isovalue and press the Submit button on the bottom right of the window. The isosurface will now appear on the main window.
Using an isovalue one order of magnitude smaller than the maximum value, the obtained plot should be similar to the following one:
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. Here, we modify inside the existing file input.xml the second index in the element pointstatepair as follows:
...
<kstlist>
<pointstatepair>1 7</pointstatepair>
</kstlist>
...
%%bash
mkdir -p CO2/HOMO_2
Now we run again exciting using the command
%%bash
cd CO2/HOMO_2
python3 -m excitingscripts.execute.single
cd ../../
and, after the run is completed, move to the directory 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.
%%bash
cd CO2/HOMO_2
python3 -m excitingscripts.convert_xml2xsf -d 3D -f WF3D.xml
python3 -m excitingscripts.convert_xml2xsf -d 3D -f input.xml
cat input.xsf WF3D.xsf > HOMO_2.xsf
xcrysden --xsf HOMO_2.xsf >/dev/null 2>&1 &
cd ../../
The obtained plot should be similar to the following one:
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.
Please notice that, , due to the degeneracy of the HOMO_1 and HOMO_2 orbitals, their plots could appear different from what is shown here. However, the difference will consist in the application of an artificial rotation of both orbitals around the axis of the molecule. Of course, this difference has no effect on the physical properties of the molecule.
iii) Visualization of the LUMO
Finally, we also calculate and plot the lowest unoccupied MO (LUMO) of CO2. The procedure is now known. 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:
%%bash
mkdir -p CO2/LUMO
%%bash
cd CO2/LUMO
python3 -m excitingscripts.execute.single
cd ../../
Repeating the standard visualization procedure inside the directory LUMO, we obtain the file LUMO.xsf which originates a plot which should be similar to the following image (using an isovalue around 0.015).
%%bash
cd CO2/LUMO
python3 -m excitingscripts.convert_xml2xsf -d 3D -f WF3D.xml
python3 -m excitingscripts.convert_xml2xsf -d 3D -f input.xml
cat input.xsf WF3D.xsf > LUMO.xsf
xcrysden --xsf LUMO.xsf >/dev/null 2>&1 &
cd ../../
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?
Visualize the HOMO and LUMO for the water and pyridine molecule, too.