Tutorial:Nitrogen atom

From OctopusWiki
Jump to: navigation, search

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.

The input files

This sample input file lets us obtain the ground state of the nitrogen atom, within the LDA approximation, in a closed-shell (unpolarized) configuration (as explained below, you need an auxiliary .xyz input). 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
UnitsOutput = eV_Angstrom

Nitrogen_mass = 14.0

 'N' | species_pseudo | db_file | 'PSF/N.psf' | lmax | 1 | lloc | 0 | mass | Nitrogen_mass

XYZCoordinates = 'N.xyz'

ExtraStates = 1
  2 | 1 | 1 | 1

BoxShape = sphere
Radius = 5.0*angstrom
Spacing = 0.18*angstrom

The most important new variables are:

  • UnitsOutput = eV_Angstrom: Two different unit systems may be used for output: the usual atomic units (which is the default, and the ones used internally in the code); and the system in which the the Ångström is substituted for the atomic unit of length, and the electronvolt is substituted for the atomic unit of energy. You can find a more detailed description of units in Octopus in the Units page of the manual.
  • The following entry in the input file are not a variable that Octopus will read directly, but rather illustrate the possibility of writing "user-defined" values and expressions to simplify the input file. In this case, we define the nitrogen mass (Nitrogen_mass = 14.0) (note that in this case, as an exception, the value is expected to be in the so-called "atomic mass units", rather than in "atomic units"). This definition may be used elsewhere in the input file.
  • The Species block should contain the list of species that are present in the system to be studied, . In this case we have only one species: nitrogen. The first field is a string that defines the name of the species, "N" in this case. The second field defines the type of species, in this case species_pseudo . Then a list of parameters follows. The parameters are specified by a first field with the parameter name and the field that follows with the value of the parameter. Some parameters are specific to a certain species while others are accepted by all species. In our example db_file instructs Octopus to find a N.psf (a Troullier-Martins pseudopotential (SIESTA format) for this nitrogen atom) in the directory share/octopus/pseudopotentials/PSF. Then come maximum lmax - component of the pseudopotential to consider in the calculation, and the lloc - component to consider as local. Generally, you want to set the maximum l to the highest available in the pseudopotential and the local l equal to the maximum l. Finally, the mass of the species can also be modified from the default values by setting mass parameter.
  • 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:
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 in each 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, which we get with this line in the input file.
  • %Occupations block: Generally, 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 in this 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).


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

****************************** Species *******************************
Reading pseudopotential from file:
      Calculating atomic pseudo-eigenfunctions for species N ....
Info: l =  0 component used as local potential.
Info: l =  1 is maximum angular momentum considered. 
Number of orbitals: total =     16, bound =     16

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
  Octopus will run in 3 dimension(s).
  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 =     127183
 Grid Cutoff [eV] =  1160.586810    Grid Cutoff [Ry] =    85.301565

This step is about the construction of the mesh. As requested in the input file, a sphere of radius 5 Å is used, which contains a cubic regular real-space grid with spacing 0.18 Å. 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 density cutoff of 1160.595 eV = 42.6 Hartree (except that in this case there is no artificial periodic repetition of the system).

Let us now take a look at how the code pursues its calculation. After some output you should see something like:

Info: Performing initial LCAO calculation with      4 orbitals.
Info: Getting Hamiltonian matrix elements.
ETA: .......1......2.......3......4......5.......6......7.......8......9......0 

Eigenvalues [eV]
 #st  Spin   Eigenvalue      Occupation
   1   --   -17.403388       2.000000
   2   --    -6.420082       1.000000
   3   --    -6.420082       1.000000
   4   --    -6.420082       1.000000

This is the first step of a ground-state calculation: obtaining 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 written to standard output.

*********************** SCF CYCLE ITER #    1 ************************
etot = -2.54498867E+02 abs_ev   =  4.33E-01 rel_ev   =  7.95E-03
                       abs_dens =  2.82E-01 rel_dens =  5.64E-02 
Matrix vector products:    108
Converged eigenvectors:      0 

#  State  Eigenvalue [eV]  Occupation    Error
      1      -17.470027    2.000000   (1.1E-04)
      2       -6.520121    1.000000   (1.1E-04)
      3       -6.520121    1.000000   (1.1E-04)
      4       -6.520121    1.000000   (1.1E-04)

Density of states:


Elapsed time for SCF step     1:          2.51

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 rel_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 ConvRelEv).
  • 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 wavefunctions. 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.


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. You will notice that octopus did not give any warning after the line

Info: Loading restart information.

This is useful if you change slightly the parameters of the simulation (for example the XC functional or the convergence criteria). If you change the grid parameters Octopus will not be able to restart from the previous calculation. If you do not want Octopus to try to restart a calculation, you can set the variable FromScratch.

Finding a good spacing

Convergence with spacing of N

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 BoxShape = parallelepiped 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, the simplest one being to change the input file by hand and run Octopus each time. But we can use this little bash script. Make a file called spacing.sh and then type source spacing.sh (NOT sh spacing.sh or bash spacing.sh!):

echo "#Sp    Energy        s_eigen   p_eigen" > spacing.log
list="0.26 0.24 0.22 0.20 0.18 0.16 0.14"
export OCT_PARSE_ENV=1 
for Spacing in $list 
 export OCT_Spacing=$(echo $Spacing*1.8897261328856432 | bc)
 octopus >& out-$Spacing
 energy=`grep Total static/info  | head -1 | cut -d "=" -f 2`
 seigen=`grep  "1   --" static/info | head -1 |  awk '{print $3}'`
 peigen=`grep  "2   --" static/info | head -1 |  awk '{print $3}'`
 echo $Spacing $energy $seigen $peigen >> spacing.log
 rm -rf restart
unset OCT_Spacing

It uses the feature that one can override variables in the inp file by defining them as environment variables in the shell. Note that we unset the variable OCT_Spacing at the end of the script, to avoid it to remaining defined in case you run the script directly in the shell. Plot your results in your favorite plotting program, e.g. gnuplot.

The results, for this particular example, are shown in the figure. In reading it, note that the most accurate results are at the left (smallest spacing). 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.

Previous Tutorial:Hydrogen atom - Next Tutorial:Methane molecule

Back to Tutorial