Tutorial:1D Helium

From OctopusWiki

1 Introduction
2 Input
3 Running
4 Exercises

Introduction

The next example will be the helium atom in one dimension which also has two electrons just as we used for the harmonic oscillator. Instead of describing two electrons in one dimension we will describe one electron in two dimensions. The calculation in this case is not a DFT one, but an exact solution of the Schrödinger equation -- not an exact solution of the helium atom, however, since it is a one-dimensional model. Calling $x\,$ and $y\,$ the coordinates of the two electrons, the Hamiltonian would be (note that the usual Coulomb interaction between particles is usually substituted, in 1D models, by the soft Coulomb potential, $u(x)=(1+x^2)^{(-1/2)}\,$ ):

$\hat{H} = -\frac{1}{2}\frac{\partial^2}{\partial x^2} -\frac{1}{2}\frac{\partial^2}{\partial y^2} +\frac{-2}{\sqrt{1+x^2}}+\frac{-2}{\sqrt{1+y^2}}+\frac{1}{\sqrt{1+(x-y)^2}}.$

Instead of describing two electrons in one dimension, however, we may very well think of one electron in two dimensions, subject to a external potential with precisely the shape given by:

$\frac{-2}{\sqrt{1+x^2}}+\frac{-2}{\sqrt{1+y^2}}+\frac{1}{\sqrt{1+(x-y)^2}} \,,$

Since the Hamiltonian is identical, we will get the same result. Whether we regard $x\,$ and $y\,$ as the coordinates of two different particles in one dimension or as the coordinates of the same particle along the two axes in two dimensions is entirely up to us. Since it is usually easier to treat only one particle, we will solve the one-dimensional Helium atom in two dimensions. We will also therefore get a "two-dimensional wave-function". In order to plot this wave-function we specify an output plane instead of an axis.

[edit]

Input

With the different potential and one more dimension the new input file looks like the following

CalculationMode = gs

Dimensions = 2
TheoryLevel = independent_particles

BoxShape = parallelepiped
lsize = 8
spacing = 0.1

Output = wfs
OutputHow = plane_z
#OutputWfsNumber = "1-2,4,6"

%Species
  "helium" | 4 | spec_user_defined | 1 | "-2/(1+x^2)^(1/2)-2/(1+y^2)^(1/2)+1/(1+(x-y)^2)^(1/2)"
%

%Coordinates
  "helium"| 0 | 0 | 0
%

We have also commented out the line that chooses the wave-functions for the output. (Octopus will always ignore anything that comes after a #.) Since we went back to CalculationMode = gs only the ground-state is calculated.

We named the species "helium" instead of "He" because "He" is already the name of a pseudopotential for a 3-d helium atom.

[edit]

Running

The calculation should converge within 8 iterations. If it happens that you want to stop the calculation and then restart it there is good news, Octopus allows you to do this. If you now just rerun Octopus the calculation will converge with just one iteration, because you use the already converged data from restart/gs. Results are summarized in the file static/info, where you can find

...
Independent Particles
SCF converged in    8 iterations

Eigenvalues [H]
 #st  Spin   Eigenvalue     Occupation
   1   --    -2.238257       1.000000

Energy:
      Total       =     -2.23825730
      Ion-ion     =      0.00000000
      Eigenvalues =     -2.23825730
      Hartree     =      0.00000000
      Int[n v_xc] =      0.00000000
      Exchange    =      0.00000000
      Correlation =      0.00000000
      Kinetic     =      0.28972533
      External    =     -2.52798263

As we are running with non-interacting electrons, the Hartree, exchange and correlation components of the Total energy are zero. Also the ion-ion term is zero, as we only have one "ion".

Now you can do just the same thing we did for the harmonic oscillator and change the CalculationMode to unocc and calculate five unoccupied states.

The results of this calculation can be found in the file static/eigenvalues. In this case it looks like

Ground-state of He in 1D

1st excited-state of He in 1D

2nd excited-state of He in 1D

Some of the unoccupied states are not fully converged!
Criterion =      0.100000E-05

Eigenvalues [H]
 #st  Spin   Eigenvalue     Occupation       Error
   1   --    -2.238257       1.000000      (0.0E+00)
   2   --    -1.815718       0.000000      (9.4E-07)
   3   --    -1.701549       0.000000      (2.1E-05)
   4   --    -1.629240       0.000000      (3.1E-08)
   5   --    -1.608656       0.000000      (3.1E-08)
   6   --    -1.509599       0.000000      (2.9E-08)

Note that the 3rd state did not converge to our criterion. This is not a big problem for us, but please be attentive to this kind of warning. Furthermore, the higher-lying states are still not converged with the size of the box. You will hear about convergence in another tutorial.

Apart from the eigenvalues and occupation numbers we asked Octopus to output the wave-functions. To plot them, we will use the gnuplot, you start it and type

$ gnuplot
gnuplot> set hidden3d
gnuplot> set pm3d
gnuplot> set contour
gnuplot> set ticslevel 0
gnuplot> unset key
gnuplot> unset surface
gnuplot> splot 'static/wf-k001-st0001-sp1.z=0' using 1:2:3
gnuplot> set term x11 1
gnuplot> splot 'static/wf-k001-st0002-sp1.z=0' using 1:2:3
gnuplot> set term x11 2
gnuplot> splot 'static/wf-k001-st0003-sp1.z=0' using 1:2:3

The first 6 lines are for eye-candy purposes. Then we plot the ground-state, 1st and 2nd excited-state wave-functions. Which correspond to singlet and which to triplet states?

[edit]

Exercises

Calculate the helium atom in 1D assuming that the 2 electrons of helium do not interact (using Dimensions = 1). Can you justify the differences?

See how the results change when you change the interaction. Often one models the Coulomb interaction by $1/\sqrt{a^2+r^2}\,$ , and fits the parameter $a\,$ to reproduce some required property.