Tutorial:1D Harmonic Oscillator
From OctopusWiki
Contents |
Input
As the first example we use the standard textbook harmonic oscillator in one dimension and fill it with two non-interacting electrons. Now we just have to tell octopus to do exactly that. The program expects you to write an input file called inp, where all the initial information is stored. octopus will always look for this file in the directory you are working in. So we suggest you just make a subdirectory for every example in this tutorial and then you can always give the input file the same name.
octopus already knows a couple of things that make the input a lot easier. Most of the usual functions like sin and cos but also sqrt and many more are known. octopus also knows the coordinates x, y, z, r that you will need to put in potentials. It is also capable of working with pseudopotentials, as we will see later on. Some parts of the input are expected to be given as strings and have to be written in quotation marks.
So now it's time to write the first input file. So get your favourite text editor, write the following lines and save the file as inp.
CalculationMode= gsDimensions= 1TheoryLevel= independent_particlesradius= 10spacing= 0.1 %Species"HO" | 0.0 | spec_user_defined | 2 | "0.5*x^2" % %Coordinates"HO" | 0 | 0 | 0 %
-
CalculationMode= gs: This variable defines the run mode - please consult the manual for the full list of the possible run modes. In this case we set it to gs, which instructs the code to start a ground-state calculation. Note that gs is simply a mnemonic for the value 1, which is the one that determines the run mode. I.e. if you writeCalculationMode= 1, you get the same behavior. However, we strongly discourage the use of the numerical values, as these may change - the mnemonics are more stable. Most variables whose list of accepted values are numeric have this kind of mnemonic helps - the full list of equivalence is given in the file SHARE/octopus/variables. (Octopus is installed, as described in the manual under some directory which hereafter we will call PREFIX. The executables are then placed in PREFIX}/bin, and a directory called PREFIX/share/octopus will be created. We will call SHARE = PREFIX/share.)
-
Dimensions= 1: This means that the code should run for 1-dimensional systems. Other options are 2, or 3 (the default). At the moment, we do not have the possibility of running in 4 or higher dimensional spaces!
-
TheoryLevel= independent_particles: We tell Octopus to treat the electrons as non-interacting.
-
radius= 10.0: The radius of the sphere, obviously. As the default unit system is atomic units, this means 5 bohr.
-
spacing= 0.1: As you should know, Octopus works in a real space regular cubic mesh. This variable defines the spacing between points, a key numerical parameter, in some ways equivalent to the energy cutoff in plane wave calculations.
- %
Species: The entry is not just the definition of a variable, but rather of a full set of them - a "block" of variables. The beginning of a block is marked by the %identifier line, and ended by a % line. The identifier in this case isSpecies. This block should contain the list of species that are present in the system to be studied. So we defineSpecies, which is nothing other than the atoms or elements of the system we want to calculate. In our case it is a harmonic oscillator, so we give the default name "HO", the second entry corresponds to its mass that we choose to be zero, and the third entry indicates that this is a user-defined species (the other choices correspond to pseudopotentials). The next entry is the number of valence electrons, in our case 2. Finally, the last entry is the potential that resembles this kind of atom, so just the harmonic oscillator potential. Note that the potential and the name have to be given as strings.
- %
Coordinates: In this block we can list the "atoms" in our calculation, one per line. In this case we do not really have an atom, but a Harmonic Oscillator, but the difference is irrelevant for octopus. What the line means is to add the external potential defined for species "HO" in the position (0,0,0).
Output
Now one can execute this file by running octopus. This is what octopus outputs to stdout:
% ../trunk/src/main/octopus
<><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><>
___
.-' `'.
/ \
| ;
| | ___.--,
_.._ |0) ~ (0) | _.---'`__.-( (_.
__.--'`_.. '.__.\ '--. \_.-' ,.--'` `""`
( ,.--'` ',__ /./; ;, '.__.'` __
_`) ) .---.__.' / | |\ \__..--"" """--.,_
`---' .'.''-._.-'`_./ /\ '. \ _.-~~~````~~~-._`-.__.'
| | .' _.-' | | \ \ '. `~---`
\ \/ .' \ \ '. '-._)
\/ / \ \ `=.__`~-.
jgs / /\ `) ) / / `"".`\
, _.-'.'\ \ / / ( ( / /
`--~` ) ) .-'.' '.'. | (
(/` ( (` ) ) '-;
` '-; (-'
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2, or (at your option)
any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program; if not, write to the Free Software
Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
<><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><>
At the beginning, it just prints a pretty picture and the short version of the license.
Running octopus
Version : 3.2.0
Revision : 6135
Build time : Fri Nov 27 16:16:18 CET 2009
Configuration options : max-dim=3
Optional libraries :
Architecture : x86_64
C compiler : gcc-4.4
C compiler flags : -pg -g -Wall -pedantic -march=native -funroll-loops -O3 -ffast-math -Wstrict-aliasing=2 -fgcse-lm -fgcse-sm -ffast-math -ftree-vectorize -pipe -I/opt/gfortran//include
Fortran compiler : gfortran-4.4
Fortran compiler flags : -pg -g -Wall -pedantic -march=native -funroll-loops -O3 -ffast-math -Wstrict-aliasing=2 -fgcse-lm -fgcse-sm -ffast-math -ftree-vectorize -pipe -fexternal-blas -ffree-line-length-none -fmax-identifier-length=31 -fbounds-check -finit
The octopus is swimming in lascar (Linux)
Calculation started on 2009/11/27 at 16:18:48
Then comes some generic (self-explanatory) information. Note that it also gives you the revision number, the compiler, and the compiler flags used. You should always include this number when submitting a bug report!
************************** Calculation Mode **************************
Input: [CalculationMode = gs]
**********************************************************************
This is just telling us that the input variable CalculationMode has the value gs. In this case, this is just what you typed in the input file, but in other cases it can reflect the default value of some variable. You will find several of these lines in the output - you should check that they have the correct values.
****************************** Species *******************************
Specie "HO" is an user-defined potential.
Potential = 0.5*x^2
**********************************************************************
This is the information relative to the species. In this case there is only one called "HO" that is defined by the potential "0.5*x^2".
Input: [SpinComponents= unpolarized] Input: [CurvMethod= curv_uniform] Input: [DerivativesSpace= real_space] Input: [DerivativesStencil= stencil_star]
This section of the output contains fairly technical information that you can ignore safely for now. Just retain that the calculation is done for spin-compensated systems (spin-unpolarized).
************************** Parallelization *************************** Octopus will run in *serial* **********************************************************************
This section of the output tells us that octopus is running in serial, i.e. in a single processor. You may also run in parallel (the code must be compiled using MPI for that).
******************************** Grid ******************************** Simulation Box: Type = sphere Radius [b] = 10.000 The octopus will run in 1 dimension(s). The octopus will treat the system as periodic in 0 dimension(s). Main mesh: Spacing [b] = ( 0.100,-1.000,-1.000) volume/point [b^3] = 0.10000 # inner mesh = 201 # total mesh = 209 Grid Cutoff [H] = 493.480 **********************************************************************
This is information about the mesh used in the calculations. As you probably know, octopus is a real-space, grid-based code, which means that all quantities are discretized in real space. For now just notice that we have 201 in the "inner mesh" (the simulation box that contains our system) and 209 in the "total mesh" (which contains the inner mesh plus some points that are used to handle the boundary conditions). The total number of points depends, obviously, on the spacing between the points and on the size of the simulation box.
****************************** Hartree ******************************* Input: [PoissonSolver= direct1D] ********************************************************************** Info: Treating the electrons as non-interacting Input: [RelativisticCorrection= non_relativistic] Input: [TDGauge= length] Input: [AbsorbingBoundaries= no_absorbing]
Again a technical section. For now just retain that electrons are treated as "non-interacting", which means that the Hartree and exchange-correlation terms are not included.
******************** Loading restart information ********************* ** Warning: ** Could not load wave-functions from 'restart/gs' ** Starting from scratch!
By default, octopus tries to restart all calculations. For that it needs the wave-functions stores in the restart directory. In this case it could not find these wave-functions, so it starts from scratch (see next section)
Info: Unnormalized total charge = 2.000000
Info: Renormalized total charge = 2.000000
Info: Setting up Hamiltonian.
Input: [LCAOStart = lcao_none]
Starting from scratch means that octopus generates a starting density from the sum of atomic densities. This is then renormalized to integrate to the total number of electrons present in the system. Finally, it obtains a first estimate for the wave-functions using a linear combination of atomic orbitals (LCAO) technique. For that it needs the atomic wave-functions that are unfortunately not available for user defined potentials, so this step is skipped in this example. The starting orbitals will then be random.
Info: SCF using real wavefunctions. Input: [What2Mix= density] (what to mix during SCF cycles) Input: [TypeOfMixing= broyden] Input: [EigenSolver= cg] Input: [Preconditioner= pre_filter]
Very often one can work with real wave-functions. This is particularly helpful as calculations with real wave-functions are much faster than with complex ones. However, if a magnetic field is present, if the system is periodic, or if spin-orbit coupling is present, complex wave-functions are mandatory. But don't worry: the program is able to figure out by itself what to use.
Furthermore, during the self-consistent procedure one has to use a mixing scheme to help convergence. One can mix either the density or the potential, and there are several mixing schemes available. The eigen solver used will be simple conjugate gradients (cg), and a preconditioner is used to speed up its convergence.
Now start the self-consistent cycles. Note that, as the electrons are non-interacting, there is no self-consistency needed!
*********************** SCF CYCLE ITER # 1 ************************
etot = 1.00000000E+00 abs_ev = 3.26E-01 rel_ev = 3.26E-01
abs_dens = 3.27E+00 rel_dens = 1.64E+00
Matrix vector products: 27
Converged eigenvectors: 0
Eigenvalues [H]
#st Spin Eigenvalue Occupation Error
1 -- 0.500000 2.000000 (4.5E-05)
Elapsed time for SCF step: 0.01
**********************************************************************
*********************** SCF CYCLE ITER # 2 ************************
etot = 1.00000000E+00 abs_ev = 3.60E-11 rel_ev = 3.60E-11
abs_dens = 2.29E+00 rel_dens = 1.14E+00
Matrix vector products: 24
Converged eigenvectors: 1
Eigenvalues [H]
#st Spin Eigenvalue Occupation Error
1 -- 0.500000 2.000000 (9.0E-07)
Elapsed time for SCF step: 0.00
**********************************************************************
*********************** SCF CYCLE ITER # 3 ************************
etot = 1.00000000E+00 abs_ev = 7.66E-15 rel_ev = 7.66E-15
abs_dens = 1.44E-05 rel_dens = 7.18E-06
Matrix vector products: 2
Converged eigenvectors: 1
Eigenvalues [H]
#st Spin Eigenvalue Occupation Error
1 -- 0.500000 2.000000 (4.1E-07)
Elapsed time for SCF step: 0.00
**********************************************************************
*********************** SCF CYCLE ITER # 4 ************************
etot = 1.00000000E+00 abs_ev = 1.11E-16 rel_ev = 1.11E-16
abs_dens = 2.02E-08 rel_dens = 1.01E-08
Matrix vector products: 2
Converged eigenvectors: 1
Eigenvalues [H]
#st Spin Eigenvalue Occupation Error
1 -- 0.500000 2.000000 (5.5E-07)
Elapsed time for SCF step: 0.00
**********************************************************************
Info: SCF converged in 4 iterations
Several pieces of information are output per self-consistency step. The first line gives the total energy (etot), and the absolute (abs_ev) and relative (rel_ev) convergence in the eigenvalues. The second line gives the absolute convergence. Then come the number of Hamiltonian - wave-function products used, followed by the number of converged eigenvectors. Note that sometimes octopus says that the SCF cycle converged but the number of "Converged eigenvectors" is smaller than the total number of eigenvectors. Do not worry: your calculation is very likely converged, so you can probably ignore this.
Calculation ended on 2009/11/27 at 16:18:48
And it is over...
Files
After finishing the calculation you will find a series of files in the directory you ran:
% ls exec inp restart static
exec
This directory contains files regarding the execution of octopus:
% ls exec nl_operator_prof oct-status-finished out.oct
The first file, nl_operator_prof contains some profiling information on your system. The second indicates that Octopus has finished. Finally, out.oct is a plain-text file that contains the information octopus reads from your input file. If you look at it, you will find the information you included in your file and additional entries which have a #default comment to it. So things you did not specify in your input file will be assigned their default values. You can also use this file to find out what options you can specify and what the appropriate variable is called. So it's a good idea to have a look at this file from time to time.
restart
This is where the files needed to restart a calculation are stored. It may contain several sub-directories depending on the calculations previously performed. In this case, it just contains one:
% ls restart gs % ls restart/gs 0000000001.obf Lxyz mesh occs states wfns
octopus stores each individual state in a different binary (yet platform independent) file. In this case, we only have one state that is in the file 0000000001.obf. The other files are text files that contain diverse control information. It is unlikely that you will ever have to work directly with these files, but you may take a look around if you are curious.
static
This directory contains the results from a static (in this case ground-state) calculation.
% ls static info momentum
The file info is a plain text file so just have a look at it. It will give you detailed information about the mesh used, the results with errors and convergence criteria, and so on. This is clearly the most important file from this run!
Exercises
Now we can play a little bit with the input file and add some other features. For example we could think of not only calculating the ground state but also some unoccupied states. This can be done by changing the CalculationMode to
CalculationMode = unocc
in the input file. In this mode octopus will not only give you the ground state but also five excited states. This is just a default value that you can change yourself by putting an extra line
NumberUnoccStates = 10
that will calculate 10 excited states. A thing to note here is that octopus will need the density for this calculation. (Actually for non-interacting electrons the density is not needed for anything, but since octopus is designed for interacting electrons it will try to read the density anyways.) So if you have already performed a static calculation (CalculationMode = gs) it will just use this result.
Compared to the ground-state calculation we also have to change the convergence criterion. The unoccupied states do not contribute to the density so they might not (and actually will not) converge properly if we use the density for the convergence check. Therefore, octopus now checks whether all calculated states converge seperately by just looking at the biggest error in the wave-functions. From the calculation you get another file in the static directory called eigenvalues. The info file will only contain the information about the ground-state, all eigenvalues and occupation numbers will be in the eigenvalues file.
If we also want to plot, say, the wave-function, at the end of the calculation, we have to tell octopus to give us this wave-function and how it should do this. We just include
Output= wfsOutputHow= axis_x
The first line tells octopus to give out the wave-functions and the second line says it should do so along the x-axis. We can also select the wave-functions we would like as output, for example the first and second, the fourth and the sixth. (If you don't specify anything octopus will give them all.)
OutputWfsNumber = "1-2,4,6"
octopus will store the wave-functions in the same folder static where the info file is, under a meaningful name. They are stored as pairs of the x coordinate and the value of the wave-function at that position x. One can easily plot them with gnuplot or a similar program.
It is also possible to extract a couple of other things from octopus like the density or the Kohn-Sham potential. We will get to some of them in the following examples.
