e-H scattering

In this tutorial, we will show how to model the problem of an electron wavepacket scattering on an hydrogen atom using Octopus. In order to speed up the calculations and make it easier to plot the different quantities of interest, we will do the calculation in 1D. It should be straightforward to change this example to simulate a “real” 3D hydrogen atom.

Ground state

The first thing to do is to tell Octopus what we want it to do. We will start with an input file containing only the description of a 1D hydrogen atom, modeled by a soft Coulomb potential:


 FromScratch = yes
 CalculationMode = gs
 
 Dimensions = 1 
  
 Radius = 50
 Spacing = 0.15
  
  %Species
  "H" | species_soft_coulomb | softening | 1 | valence | 1
  %
  
  %Coordinates
  "H" | -10
  %
  
  SpinComponents = polarized
  
  ExtraStates = 1̣
  
  %Occupations
   1 | 0
   0 | 0
  %

These input variables should already be familiar as they have been explained in previous tutorials. Here we simply note the following:

Output

Now one can execute this file by running Octopus. Here are a few things worthy of a closer look in the standard output. We are doing an LDA calculation, as we want to investigate the role of electron-electron interaction here, and, as a possible application, benchmark the performance of exchange-correlation functionals. As we are doing a one dimensional calculation, we see that Octopus has selected the corresponding 1D version of the LDA functional:



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

Exchange-correlation:
  Exchange
    Exchange in 1D (LDA)
  Correlation
    Casula, Sorella & Senatore (LDA)

e-H scattering

We are now ready to perform a simulation of the e-H scattering problem. For this, we use the following input file:


 FromScratch = yes
 CalculationMode = td
 
 Dimensions = 1 
 
 Radius = 50
 Spacing = 0.15
  
 %Species
 "H" | species_soft_coulomb | softening | 1 | valence | 1
 %
  
 %Coordinates
 "H" | -10
 %
  
 SpinComponents = polarized
 
 ExtraStates = 1̣
 
 ExcessCharge = -1
 %Occupations
  1 | 1
  0 | 0
 %
 
 a = 0.1
 x0 = 10.0
 p = -1.5
 %UserDefinedStates 
   1 | 2 | 1 | formula | "exp(-a*(x-x0)^2+i*p*(x-x0))"
 %
 
 %Output
  density 
  potential
 %
 OutputFormat = axis_x
 
 TDPropagationTime = 0.8 * femtosecond

In order to introduce the electron wavepacket to scatter on the hydrogen atom, we have done the following change:

Fig. 1. e-H scattering computed at the ALDA level.
Fig. 1. e-H scattering computed at the ALDA level.

The result of the calculation is shown in Fig. 1. What we are interested in is the density of the up channel, where both electrons are, as well as the corresponding potential. These quantities are found in the file ‘output_iter/td.XXXXXXX/density-sp1.y=0,z=0’ and ‘output_iter/td.XXXXXXX/vxc-sp1.y=0,z=0’. We note a few interesting details. First, we observe no significant reflection of the wavepacket density on the hydrogen atom. This is a known deficiency of the adiabatic approximation. Moreover, towards the end of the simulation, both the density and the exchange-correlation potential show strong oscillations on the left-hand side of the simulation. This is due to the artificial reflection of the electron wavepacket at the border of the simulation box. This can be easily fixed by double the size of the box (Radius = 100), and/or by using absorbing boundaries.

Gnuplot script

In order to produce the above movie, one can use the following gnuplot script



reset 

set terminal pngcairo size 600,700


do for [t=0:1112:50] {

set output sprintf('figs/potential_%04i.png',t)

set multiplot

set xrange [-50:50]

set tmargin 1.5
set bmargin 1.5

set size 1.0, 0.6
set origin 0.0, 0.4 
set ylabel "n(x)" offset 2.0,0


au2fs = 0.024189
dt = 0.02973;

set label sprintf('t=%.2f fs', t*dt*au2fs) at -15,0.45
set yrange [0:0.5]
plot sprintf('output_iter/td.%07i/density-sp1.y=0,z=0', t) u 1:2 w l lw 2 notitle
  
unset label
set size 1.0, 0.4
set origin 0.0, 0.0 
set yrange [-0.5:0.0]
set ylabel "v_{xc}(x)" offset 2.0,0
set xlabel "x [Bohr]" offset 0, 0.6
plot sprintf('output_iter/td.%07i/vxc-sp1.y=0,z=0', t) u 1:2 w l ls 2 lw 2 notitle

unset label
unset multiplot
}

and convert it to a movie using convert

  convert -delay 20 -loop 0 figs/*.png scattering.gif