Tutorial:Nitrogen atom

From OctopusWiki

Now we will move to a more complicated (and realistic) input file. We will obtain the ground state of an atomic system. We will give a rather detailed description of the input and output files for this example.

Contents 1 The input file 2 Output 3 Restarting and Visualization 4 Finding a good spacing

The input file

The sample input file permits to obtain the ground state of the nitrogen atom, within the LDA approximation, in a closed-shell (unpolarized) configuration. Note that this is not the correct ground state of the nitrogen atom! However, it will permit us to describe some of the most important input variables:

CalculationMode = gs
units = eV_Angstrom

Nitrogen_mass = 14.0
Nitrogen_z = 7

%Species
  'N' | Nitrogen_mass | spec_ps_psf | Nitrogen_z | 1 | 0
%

XYZCoordinates = 'N.xyz'

ExtraStates = 1
%Occupations
  2 | 1 | 1 | 1
%

BoxShape = sphere
radius = 5.0
spacing = 0.18

The most important new variables are:

• units = eV_Angstrom: Two different unit systems may be used both for input and output: the usual atomic units (which is the default, and the ones used internally in the code); and the system in which the atomic unit of length is substituted by the Ångström, and the atomic unit of energy is substituted by the electronvolt. This is the case in this input file.

• The following two entries in the input file are not Octopus variables, but rather illustrate the possibility of writing "user defined" values and expressions. In this case, we define the atomic number of Nitrogen (Nitrogen_mass = 14.0), and its mass (Nitrogen_z = 7) (note that in this case, as an exception, the value is expected to be in the so called "atomic units of mass", rather than in "atomic units" or in the other 'eVA' system). These definitions may be used elsewhere in the input file.

• The Species block is a bit different from the the previous tutorials. The third entry is now spec_ps_psf. It instructs Octopus to use a Troullier-Martins pseudopotential (siesta format) for this Nitrogen atom. In practice, this means that Octopus will try to find a N.psf file either in the working directory (from where Octopus was launched) or in the directory share/octopus/PP/PSF. The following two entries require some knowledge of what nonlocal pseudopotentials are, and how Troullier-Martins pseudopotentials are constructed. We do not describe them here; it suffices to say that the two values entered here are perfectly reasonable for Nitrogen.

• XYZCoordinates = 'N.xyz': The geometry of the molecule (in this case, a single atom in the grid origin) is described in this case in a file with the well known XYZ format. The file for this outrageously simple case is given by:

1
This is a comment line
N 0 0 0

• ExtraStates = 1: By default, octopus performs spin-unpolarized calculations (restricted closed-shell, in Hartree-Fock terminology). It then places two electrons per orbital. The number of orbitals, or Kohn-Sham states, is then calculated by counting the number of valence electrons present in the system, and dividing by two. In this case, since we have five valence electrons, the code would use three orbitals. However, we know beforehand that the HOMO orbital has a three-fold degeneracy, and as a consequence we need to put each one of the three p electrons in a different orbital. We therefore need one more orbital, and this is the request of this entry in the input file.

• %Occupations block: In principle, the occupations of the Kohn-Sham orbitals are automatically decided by the code, filling the lowest-energy orbitals. However, if we have degeneracies in the LUMO as it is the case, the user may want to accommodate the electrons in a certain predefined way. In this example, the obvious way to fill the orbitals of the Nitrogen atom is to put two electrons in the first and deepest orbital (the s orbital), and then one electron on each of the second, third and fourth orbitals (the p orbitals, which should be degenerate).

• BoxShape = sphere: This is the choice of the shape of the simulation box, which in this case is set to be a sphere (other possible choices are minimum, cylinder, or parallelepiped).

Output

Once you have constructed the input file, you may launch Octopus on it. The new sections of the output are

****************************** Species *******************************
Reading pseudopotential from file:
      '/home/marques/software/octopus-svn/share/octopus/PP/PSF/N.psf'
      Calculating atomic pseudo-eigenfunctions for specie N ....
      Done.
Info: l =  0 component used as local potential
**********************************************************************

Here the code searches for the needed pseudopotential files, and informs the user about its success or failure. In this case, only the (default) N.psf file is required and processed.

******************************** Grid ********************************
Simulation Box:
  Type = sphere
  Radius  [A] =   5.000
The octopus will run in 3 dimension(s).
The octopus will treat the system as periodic in 0 dimension(s).
Main mesh:
  Spacing [A] = ( 0.180, 0.180, 0.180)    volume/point [A^3] =  0.00583
  # inner mesh =    89727
  # total mesh =   158807
  Grid Cutoff [eV] =  1160.595
**********************************************************************

This step regards the construction of the mesh. As requested in the input file, a sphere of radius 5 Ångström is used, which contains a cubic regular real-space grid with 0.18 Ångström of spacing. This implies 89727 points (inner mesh = 89727). For the sake of comparison with plane-wave-based codes, this is more or less equivalent to a plane-wave calculation that imposes a cutoff of 1160.595 eV = 42.6 Hartree (except that in this case there is no artificial periodic repetition of the system).

****************************** Hartree *******************************
Input: [PoissonSolver = isf]
 Build the kernel using a sum of 89 gaussians
 Do a 3D HalFFT for the kernel
**********************************************************************

The next four lines are concerned with the Poisson solver - since no input variable regarding this issue was given, the code makes use of the default option. This option is usually good enough for normal molecules.

************************ Exchange-Correlation ************************
Exchange and correlation:
  Exchange
    Slater exchange (LDA)
    [1] P.A.M. Dirac, Proceedings of the Cambridge Philosophical Society 26, 376 (1930)
    [2] F. Bloch, Zeitschrift fuer Physik 57, 545 (1929)
  Correlation
    Perdew & Zunger (Modified) (LDA)
    [1] Perdew and Zunger, Phys. Rev. B 23, 5048 (1981)
    [2] Modified to improve the matching between the low and high rs parts

Input: [SICCorrection = sic_none]
**********************************************************************

These show the (default) options for the exchange and correlation terms of the Hamiltonian: the local density approximation (LDA) for both exchange and correlation (the non-relativistic expression for the former case, and the parameterization of Perdew and Zunger for the latter). We are also not using any kind of self-interaction correction.

Let us now take a look at how the code pursues its calculation:

Input: [LCAOStart = lcao_full]
Info: Performing initial LCAO calculation with    8 orbitals.
Eigenvalues [eV]
 #st  Spin   Eigenvalue     Occupation
   1   --   -17.425434       2.000000
   2   --    -6.325377       1.000000
   3   --    -6.325377       1.000000
   4   --    -6.325377       1.000000

The first step is to obtain a reasonably good starting density and Kohn-Sham orbitals to feed in the self-consistent (SCF) procedure. For this purpose, Octopus performs an initial calculation restricted to the basis set of atomic orbitals (Linear Combination of Atomic Orbitals, LCAO). The resulting eigenvalues of this calculation are cast to standard output.

*********************** SCF CYCLE ITER #    1 ************************
 etot = -2.61943592E+02 abs_ev   =  1.66E-01 rel_ev   =  3.07E-03
                        abs_dens =  5.77E-02 rel_dens =  1.15E-02
Matrix vector products:    108
Converged eigenvectors:      0
Eigenvalues [eV]
 #st  Spin   Eigenvalue     Occupation       Error
   1   --   -17.482211       2.000000      (7.3E-05)
   2   --    -6.406929       1.000000      (5.3E-05)
   3   --    -6.406929       1.000000      (5.3E-05)
   4   --    -6.406929       1.000000      (5.3E-05)

Elapsed time for SCF step:          2.26
**********************************************************************

Now the SCF cycle starts. For every step, Octopus outputs several pieces of information:

• The values abs_dens and rel_dens are to monitor the absolute and relative convergence of the density, while rel_ev and abs_ev are two alternative measures of the convergence, based on measuring the difference between input and output eigenvalues. The SCF procedure, by default, is stopped when abs_dens is smaller than 10 ^{− 5}. This may be altered with the appropriate input variables (see, in the manual, the variables ConvAbsDens, ConvRelDens, ConvAbsEv and RelAbsEv).

• The line Matrix vector products: 108 tells us that the Hamiltonian was applied 108 times. This gives us an idea of the computational cost.

• The line Converged eigenvectors: 0 tells us that upon completion of the diagonalization procedure, none of the orbitals met the required precision criterion for the wave-functions. In a following example, we will modify this criterion in the input file.

• The list of eigenvalues is then printed, along with their errors: how much they deviate from "exact" eigenvalues of the current Hamiltonian. This number is the so-called "residue".

You can now take a look at the file static/info that will hold a summary of the calculation.

Restarting and Visualization

Any ground-state calculation may be restarted later (to refine it if it did not converge properly, or with any other purpose), provided that the contents of the restart directory are preserved. You can try this now, just by running Octopus again.

Also, please add the following three lines to the inp file:

Output = wfs
OutputWfsNumber = "1-4"
OutputHow = dx + axis_x + axis_y + axis_z

You will now notice that the convergence procedure is stopped in only one iteration; the reason is that the starting point was now the (already converged) output state of the previous run. This is why the standard output has the following lines:

******************** Loading restart information *********************
All the needed files were succesfully read.
**********************************************************************

Now take a look at the static directory. Besides the info file there are a bunch of new files, called wf-001-00x-1.dx, and wf-001-00x-1.a=0,b=0, where x runs from 1 to 4 (the four KS states), and where a and b are either x, y or z. Instructing Octopus to generate these files was the task of the input variables that you just added. These files contain various functions related to the system (in this case, wavefunctions) for their visualization.

• Output = wfs: This variable asks Octopus to print the wavefunctions. One may also wish to print densities, potentials, etc. In the manual you may find the corresponding variables (density, etc).

• OutputWfsNumber = '1-4': This variable specifies which wavefunctions to print: '1-4' asks for the KS orbitals running from one to four; '1,2,7-10' would ask for the first two, and the ones from seven to ten.

• Output = dx + axis_x + axis_y + axis_z This variable tells the code to print the functions in the format native to the OpenDX program. This program permits "sophisticated" data visualization (isosurfaces, contour plots, etc). In the figure you may see isosurfaces of the obtained p_x, p_y and p_z Kohn-Sham orbitals of Nitrogen. You will also find the files wf-001-00x.x=0,y=0, wf-001-000x.x=0,z=0 and wf-001-000x.y=0,z=0. These files contain the values of the wavefunctions along the z, y and x, respectively, in the format coordinate, real value and imaginary value. This way you can plot the functions with easier-to-use plotting programs.

Finding a good spacing

The key parameter of a real-space calculation is the spacing between the points of the mesh. In the current version of octopus, the grid is regular, so there is only one grid spacing. (In fact, if you use the parallelepiped shape for your simulation box, you may define different spacings in each direction, by using the %Spacing block, instead of the Spacing variable.) The first step in any calculation should then be making sure that this spacing is good enough for our purposes. This should be done through a convergence study, very similar to the ones performed in plane-wave calculations.

The needed spacing essentially depends on the pseudopotentials that are being used. The idea is to repeat a series of ground-state calculations, with identical input files except for the grid spacing. There are many different ways of doing it, but for now we can use this little bash script:

#!/bin/bash
list="0.26 0.24 0.22 0.20 0.18 0.16 0.14"
export OCT_PARSE_ENV=1
for OCT_Spacing in $list
do
  export OCT_Spacing
  octopus >& out-$OCT_Spacing
  energy=`grep Total static/info  | head -1 | cut -d "=" -f 2`
  echo $OCT_Spacing $energy
  rm -rf restart
done

It uses the feature that one can override variables in the inp file by defining them from the command line. The results, for this particular example, are shown in the figure. A rather good spacing for this nitrogen pseudopotential seems to be 0.18 Å. However, as we are usually not interested in total energies, but in energy differences, probably a larger one may also be used without compromising the results.