Tutorial:1D Harmonic Oscillator

From OctopusWiki
Jump to: navigation, search

Octopus has the unusual feature for DFT codes of being able to handle "model systems," i.e. ones that have a user-defined arbitrary potential as opposed to a system of real atoms. This can be useful for simplified test calculations or for computations in the "effective mass" or "envelope function" approximation, e.g. for quantum dots or quantum wells.


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. Write the following lines and save the file as inp.

FromScratch = yes
CalculationMode = gs

Dimensions = 1 
TheoryLevel = independent_particles

Radius = 10
Spacing = 0.1

  "HO" | species_user_defined | potential_formula | "0.5*x^2" | valence | 2

  "HO" | 0
  • FromScratch = yes: Do not use (or expect) restart information.
  • Dimensions = 1: This means that the code should run for 1-dimensional systems. Other options are 2, or 3 (the default). You can actually run in 4D too if you have compiled with the configure flag --max-dim=4.
  • TheoryLevel = independent_particles: We tell Octopus to treat the electrons as non-interacting.
  • Radius = 10.0: The radius of the 1D "sphere," i.e. a line; therefore domain extends from -10 to +10. As the default unit system is atomic units, this means 10 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 species name is "HO", then the potential formula is given, and finally the number of valence electrons. See Manual:Input file for a description of what kind of expressions can be given for the potential formula.
  • %Coordinates: add the external potential defined for species "HO" in the position (0,0,0).


Now one can execute this file by running Octopus. In the standard output you will find the listing of the species:

****************************** Species *******************************
Species "HO" is an user-defined potential.
   Potential = 0.5*x^2
Number of orbitals:      5

The potential is V(x) = 0.5 x^2, and 5 Hermite-polynomial orbitals are available for LCAO (the number is based on the valence). The theory level is as we requested:

**************************** Theory Level ****************************
Input: [TheoryLevel = independent_particles]

The electrons are treated as "non-interacting", which means that the Hartree and exchange-correlation terms are not included. This is usually appropriate for model systems, in particular because the standard XC approximations we use for actual electrons are not correct for "effective electrons" with a different mass.

Input: [MixField = none] (what to mix during SCF cycles)

**************************** Eigensolver *****************************
Input: [EigenSolver = cg]
Input: [Preconditioner = pre_filter]
Input: [SubspaceDiagonalization = standard]

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. Since we are using independent particles (and only one electron) we don't have to mix anything though. The eigensolver used will be simple conjugate gradients (cg), and a preconditioner is used to speed up its convergence.

Input: [LCAOStart = lcao_none]
Info: Unnormalized total charge =      2.000000
Info: Renormalized total charge =      2.000000
Info: Setting up Hamiltonian.
Orthogonalizing wavefunctions.

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. For atomic systems, the default is to find a first estimate for the wave-functions using a linear combination of atomic orbitals (LCAO) technique, using the atomic wavefunctions from the pseudopotentials. However, we do not necessarily have such corresponding wavefunctions for a user-defined potential, so LCAO is turned off by default here. The starting orbitals will then be random but orthogonal.

Info: SCF using real wavefunctions.

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.

Now we start the self-consistent cycles. Note that, as the electrons are non-interacting, there is actually no self-consistency needed.

Info: Starting SCF iteration.
ETA: .......1......2.......3......4......5.......6......7.......8......9......0

*********************** SCF CYCLE ITER #    1 ************************
 etot =  1.00000288E+00 abs_ev   =  2.26E+00 rel_ev   =  2.26E+00
                        abs_dens =  3.27E+00 rel_dens =  1.64E+00
Matrix vector products:     27
Converged eigenvectors:      0

#  State  Eigenvalue [H]  Occupation    Error
      1        0.500001    2.000000   (3.5E-03)

Elapsed time for SCF step     1:          0.00

ETA: .......1......2.......3......4......5.......6......7.......8......9......0

*********************** SCF CYCLE ITER #    2 ************************
 etot =  1.00000000E+00 abs_ev   =  2.88E-06 rel_ev   =  2.88E-06
                        abs_dens =  3.23E-03 rel_dens =  1.62E-03
Matrix vector products:     27
Converged eigenvectors:      0

#  State  Eigenvalue [H]  Occupation    Error
      1        0.500000    2.000000   (9.8E-06)

Elapsed time for SCF step     2:          0.00

ETA: .......1......2.......3......4......5.......6......7.......8......9......0

*********************** SCF CYCLE ITER #    3 ************************
 etot =  1.00000000E+00 abs_ev   =  1.27E-11 rel_ev   =  1.27E-11
                        abs_dens =  4.64E-06 rel_dens =  2.32E-06
Matrix vector products:     16
Converged eigenvectors:      1

#  State  Eigenvalue [H]  Occupation    Error
      1        0.500000    2.000000   (8.2E-07)

Elapsed time for SCF step     3:          0.00

             Info: Writing states. 2016/11/02 at 21:32:15

        Info: Finished writing states. 2016/11/02 at 21:32:15

Info: SCF converged in    3 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.

You can actually try performing the LCAO too by setting LCAOStart = lcao_states -- compare how many iterations and matrix-vector products (in total) are required now. Why? (Hint: what are Hermite polynomials?)


After finishing the calculation you will find a series of files in the directory you ran:

% ls
exec  inp  restart  static


This directory contains files regarding the execution of octopus:

% ls exec
messages            oct-status-finished parser.log

parser.log 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.


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
% ls restart/gs
0000000001.obf density        density.obf    grid           lxyz.obf       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.


This directory contains the results from a static (in this case ground-state) calculation.

% ls static

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!

Look at the eigenvalue and total energy. Compare to your expectation from the analytic solution to this problem!


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 occupied states (which contribute to the density) but also the unoccupied ones. Set the number of states with an extra line

ExtraStates = 10

that will calculate 10 empty 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 = wfs
OutputFormat = 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.

Previous Tutorial:Optical Spectra from Casida - Next Tutorial:Particle in a box

Back to Tutorial