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)
    [1] N. Helbig, J. I. Fuks, M. Casula, M. J. Verstraete, M. A. L. Marques, I. V. Tokatly, and A. Rubio, Phys. Rev. A 83, 0325
  Correlation
    Casula, Sorella & Senatore (LDA)
    [1] M. Casula, S. Sorella, and G. Senatore, Phys. Rev. B 74, 245427 (2006)

Input: [SICCorrection = sic_none]
**********************************************************************

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
%
RestartFixedOccupations = no

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:

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


set terminal pngcairo font "Monospace-Bold,20" size 600,700

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

 set output sprintf('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 potential_*.png scattering.gif
Fig. 1. e-H scattering computed at the ALDA level.
Fig. 1. e-H scattering computed at the ALDA level.