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 |
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= gsunits= 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 %Occupations2 | 1 | 1 | 1 %BoxShape= sphereradius= 5.0spacing= 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 the Ångström is substituted for the atomic unit of length, and the electronvolt is substituted for the atomic unit of energy. This is the case in this input file. You can find a more detailed description of units in Octopus in the Units page of the manual.
- The following two entries in the input file are not variables 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 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 mass units", rather than in "atomic units" or in the other 'eVA' system). These definitions may be used elsewhere in the input file.
- The
Speciesblock 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 entry gives the name of species, "N" in this case. The second one is the mass, which we conveniently defined before. The third entry is spec_ps_psf, which 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 (where Octopus was launched) or in the directory share/octopus/PP/PSF. The next field is the atomic number of the species, and then come two entries which 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 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.
- %
Occupationsblock: 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).
Output
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:
'/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 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:
Input: [LCAOStart = lcao_full]
Info: Performing initial LCAO calculation with 4 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
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.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,ConvAbsEvandRelAbsEv).
- 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.
Restarting
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 FromScratch.
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, the simplest one being to change the input file by hand and run Octopus each time. But 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 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 remain defined in case you run the script directly in the shell.
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.
