Tutorial:1D Harmonic Oscillator
As a first example we use the standard textbook harmonic oscillator in one dimension and fill it with two noninteracting electrons.
Input
The first thing to do is to tell Octopus what we want it to do. Write the following lines and save the file as inp.
FromScratch
= yesCalculationMode
= gsDimensions
= 1TheoryLevel
= independent_particlesRadius
= 10Spacing
= 0.1 %Species
"HO"  species_user_defined  potential_formula  "0.5*x^2"  valence  2 % %Coordinates
"HO"  0 %
Most of these input variables should already be familiar. Here is a more detailed explanation for some of the values:

Dimensions
= 1: This means that the code should run for 1dimensional systems. Other options are 2, or 3 (the default). You can actually run in 4D too if you have compiled with the configure flag maxdim=4.

TheoryLevel
= independent_particles: We tell Octopus to treat the electrons as noninteracting.

Radius
= 10.0: The radius of the 1D "sphere," i.e. a line; therefore domain extends from 10 to +10 bohr.
 %
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). The coordinates used in the potential formula are relative to this point.
Output
Now one can execute this file by running Octopus. Here are a few things worthy of a closer look in the standard output.
First one finds the listing of the species:
****************************** Species
*******************************
Species "HO" is an userdefined potential.
Potential = 0.5*x^2
Number of orbitals: 5
**********************************************************************
The potential is , and 5 Hermitepolynomial 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 "noninteracting", which means that the Hartree and exchangecorrelation 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)
Since we are using independent particles (and only one electron) there is no need to use a mixing scheme to accelerate the SCF 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 wavefunctions 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 userdefined potential, so LCAO is turned off by default here. The starting orbitals will then be random but orthogonal.
Now the selfconsistent cycle starts and you should see the following output:
Info: Starting SCF iteration. ETA: .......1......2.......3......4......5.......6......7.......8......9......0 *********************** SCF CYCLE ITER # 1 ************************ etot = 1.24484521E+00 abs_ev = 9.49E+01 rel_ev = 7.62E+01 ediff = 1.24E+00 abs_dens = 3.27E+00 rel_dens = 1.63E+00 Matrix vector products: 27 Converged eigenvectors: 0 # State Eigenvalue [H] Occupation Error 1 0.622423 2.000000 (6.3E+00) Elapsed time for SCF step 1: 0.00 ********************************************************************** ETA: .......1......2.......3......4......5.......6......7.......8......9......0
and after a few iterations it converges:
*********************** SCF CYCLE ITER # 22 ************************ etot = 1.00000001E+00 abs_ev = 1.95E08 rel_ev = 1.95E08 ediff = 1.95E08 abs_dens = 1.61E05 rel_dens = 8.07E06 Matrix vector products: 27 Converged eigenvectors: 0 # State Eigenvalue [H] Occupation Error 1 0.500000 2.000000 (1.4E03) Elapsed time for SCF step 22: 0.00 ********************************************************************** Info: Writing states. 2018/08/08 at 14:13:39 Info: Finished writing states. 2018/08/08 at 14:13:39 Info: SCF converged in 22 iterations
Note that, as the electrons are noninteracting, there is actually no selfconsistency needed. Nevertheless, Octopus takes several SCF iterations to converge. This is because it takes more than one SCF iteration for the eigensolver to converge the wavefunctions and eigenvalues. To improve this we can try having a better initial guess for the wavefunctions by turning the LCAO on. You can do this by setting LCAOStart
= lcao_states  compare how many iterations and matrixvector products (in total) are required now. Why? (Hint: what are Hermite polynomials?)
Exercises
The first thing to do is to look at the static/info file. Look at the eigenvalue and total energy. Compare to your expectation from the analytic solution to this problem!
We can also play a little bit with the input file and add some other features. For example, what about the other eigenstates of this system? To obtain them we can calculate 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 noninteracting 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 groundstate 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 separately by just looking at the biggest error in the wavefunctions. From the calculation you get another file in the static directory called eigenvalues. The info file will only contain the information about the groundstate; all eigenvalues and occupation numbers will be in the eigenvalues file.
If we also want to plot, say, the wavefunction, at the end of the calculation, we have to tell Octopus to give us this wavefunction and how it should do this. We just include
Output
= wfsOutputFormat
= axis_x
The first line tells Octopus to give out the wavefunctions and the second line says it should do so along the xaxis. We can also select the wavefunctions 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
= "12,4,6"
Octopus will store the wavefunctions in the same folder static where the info file is, under a meaningful name. They are stored as pairs of the xcoordinate and the value of the wavefunction 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 KohnSham potential.