
Example prepared by Javier Junquera (javier.junquera@unican.es)

WARNING: BEFORE RUNNING A CALCULATION FOR PRODUCTION,
TEST THE PSEUDOPOTENTIAL AND BASIS SETS, AND PERFORM
THE CONVERGENCE TESTS FOR THE MESH CUTOFF AND THE SAMPLING IN K-SPACE.

IN THE PRESENT EXAMPLES, AND IN ORDER TO SPEED THE CALCULATIONS,
WE PROVIDE SOME VALUES OF THESE PARAMETERS FOR YOU.

WE DO NOT WARRANTY THAT THE VALUES OF THESE PARAMETERS ARE CONVERGED.

* Edit the input file, h2o.defaultcell.fdf, and study the different variables.

* Run the code for this default input file, and dump the output into a file

siesta < h2o.defaultcell.fdf > h2o.defaultcell.fdf &
tail -f h2o.out

* Edit the output file and search for the default unit cell generated by SIESTA

You should find something like this

outcell: Unit cell vectors (Ang):
        7.286412    0.000000    0.000000
        0.000000    5.746952    0.000000
        0.000000    0.000000    5.621012

outcell: Cell vector modules (Ang)   :    7.286412    5.746952    5.621012
outcell: Cell angles (23,13,12) (deg):     90.0000     90.0000     90.0000
outcell: Cell volume (Ang**3)        :    235.3780

* Now, we will introduce our own lattice constant and lattice vectors.
Edit the file h2o.usercell.fdf and study the meaning of the new input variables:

LatticeConstant
LatticeVectors
DM.UseSaveDM

* Run siesta with the new cell introduced by the user

siesta < h2o.usercell.fdf > h2o.10Ang.out &
tail -f h2o.10Ang.out

* Change the lattice constant from 10 to 17 Ang, with steps of 1 Ang.
Save each output in a different file.

siesta < h2o.usercell.fdf > h2o.11Ang.out &
...
siesta < h2o.usercell.fdf > h2o.17Ang.out &

* Save in a file the data needed to plot a figure energy 
versus lattice constant.
In order to get this, type:

grep "Total =" h2o.*Ang.out > convergencesizecell.dat

And edit the file to delete all the non-needed information in each line.
The final file should look like this:

10    -465.798674
11    -465.798203
12    -465.797956
13    -465.797700
14    -465.797421
15    -465.797305
16    -465.797255
17    -465.797130

* Plot the figure Energy versus lattice constant

gnuplot
plot 'convergencesizecell.dat' u 1:2 w l

Can you explain why this figure has this shape?.
Would we obtain the same behaviour for a molecule of methane?. Why?.
How could you obtain a reasonable value for the energy of an
isolated water molecule with periodic boundary conditions?
