BASIS SETS I : MOLECULES


This exercise is intended to illustrate the definition of basis sets
with different numbers of orbitals and localization radii in
SIESTA. It is based on the H2O molecule. You will first use standard
basis sets generated by SIESTA, and then a basis set specifically
optimized for the H2O molecule. The ‘tips’ at
the end of this text will help you in analyzing the results of the
calculations.


Part 1: Standard basis sets

You will find three directories: SZ, DZ and DZP, in which you will
find input files for a H2O molecule, which will be solved with three
different basis sets: single-Z, double-Z and double-Z plus
polarizarion, respectively. The runs will consist in structural
relaxations via Conjugate gradients, to find the equilibrium structure
of the molecule for each of the bases utilized.

The bases are defined through these three lines:

PAO.BasisSize DZP 
PAO.EnergyShift 500.0 meV 
PAO.Splitnorm 0.15

where the first one defines the number of orbitals in the basis set
(in this case, a double-Z plus polarization, which means two shells of
s orbitals, two shells of p orbitals, and a polarization shell of d
orbitals for oxygen, and two shells of s orbitals and one shell of
polarization p orbitals for H); the second line indicates the energy
shift parameter, that determines the cutoff radius of each of the
orbitals; and the third one is the Split norm parameter, which defines
the radius of the second-Z orbital in the case of double-Z bases.

You should do runs for each of the basis sets, changing the
PAO.EnergyShift parameter, and look at the results as a function of
basis size and localization radius (or energy shift).  In particular,
you should look at:

Total energy 
Bond lenghts at the relaxed structure 
Bond angles at the relaxed structure 
CPU time 
Radius of each of the orbitals 
Shape of the orbitals

You have some files with results that you should be able to reproduce
in the Out directories.

Try to answer these questions:

What is the best basis set that you have found? Why?  How do the
results compare with experiment?  What do you consider a reasonable
basis for the molecule, if you need an accuracy in the geometry of
about 1%??  In order to assess convergence with respect to basis set
size, should you compare the results with the experimental ones, or
with those of a converged basis set calculation?


Part 2: Optimized basis sets.

You will now do the same calculations as in Part 1, but now using a
basis set specifically optimized for the H2O molecule. You can find
the files in directory DZP-opt.

Read the block which defines the basis set. Compare the parameters
with those of the calculations in Part 1.

Do the runs with this basis set, and compare the results with the
previous ones.

Tips:

Tip 1.: To find the bond lenghts, you can look in the files with 
extensions .BONDS or .BONDS_FINAL.

The bond lengths and bond angles, could also be found with the Molekel
program.  Load the h2o.ANI file (which contains the coordinates of the
molecule during the relaxation), and use the 'Geometry' option (left
button of the mouse) to measure bond lengths and angles.


Tip 2.: You can see the radii of the orbitals in the output file, just
after the line reading: printput: Basis input
-----------------------------------------------

Tip 3.: You can look at the shape of the orbitals by plotting the
contents of the .ion files produced by Siesta. These files are not
easily readable, but the 'ioncat' program can extract the relevant 
pieces of information, and the "ionplot" script can drive 'ioncat' to
plot the desired graphs. For example:

        ionplot -o 1 O

will plot the orbital with number "1" in the O.ion file.

        ioncat -i O

will print the numbers of the representative orbitals of each nlz
shell (i.e., disregarding the 'm' quantum number, which does not
affect the radial part).

        ioncat -o 1 O 

will output the data for the first orbital in O.ion.

Type "ioncat -h" to see the full set of options for the program.

The script "ionplot" accepts some of the options (such as -o)
and drives also gnuplot to plot the information.

By combining the two one could plot all the orbitals:

 for i in $(ioncat -i O); do ionplot -o $i; done

Each orbital will appear in a different window.
