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.

!mkdir -p 6-e-H-scattering

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:

%%writefile 6-e-H-scattering/inp

stdout = 'stdout_gs.txt'
stderr = 'stderr_gs.txt'

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
%

Writing 6-e-H-scattering/inp

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

Spacing: We employed a large box. This is needed for the time-dependent simulation, and, as we will see, is in fact not large enough for correctly describing our e-H problem for even 1 fs.
ExtraStates: We requested one extra state. This is not really needed for the ground-state calculation, but is crucial for the time-dependent run bellow, as there we will replace this unoccupied state by the electron wavepacket.
Occupations: We have fixed the occupations to be one in the first spin channel. This is needed for the TD run, as we will see, and also because otherwise the code might get to a different ground state than the one we are interested.

!cd 6-e-H-scattering && octopus

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 for an soft-Coulomb interaction (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, 032503 (2011)
  Correlation
    Casula, Sorella & Senatore (LDA)
    [1] M. Casula, S. Sorella, and G. Senatore.,  Phys. Rev. B 74, 245427 (2006)

**********************************************************************

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:

%%writefile 6-e-H-scattering/inp

stdout = 'stdout_td.txt'
stderr = 'stderr_td.txt'

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
TDTimeStep = 0.02973

Overwriting 6-e-H-scattering/inp

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

ExcessCharge: We introduce an extra electron in our simulation by setting this variable to -1.
Occupations: We force the extra electron to be on the second state of the first spin channel.
RestartFixedOccupations” = no: Note, that we need to override the default behaviour to keep the occupations from the restart file.
UserDefinedStates: We are replacing the originally unoccupied state obtained after the GS calculation by a user defined state, which is a Gaussian of the form \( \phi(x) = e^{-\alpha (x-x_0)^2+i p (x-x_0)}, \) where \(\alpha\) is taken to be \(\alpha=0.1\) here, \(x_0\) is the initial position of the wavepacket, taken to be at \(x_0=-10\) a.u., and \(p\) its velocity, taken to be \(p=-1.5\) a.u..

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 doubling the size of the box (Radius = 100), and/or by using absorbing boundaries.

!cd 6-e-H-scattering && octopus

from postopus import Run
import holoviews as hv
hv.extension("matplotlib")

run = Run("6-e-H-scattering")
nx = run.td.density_sp1()
vxc = run.td.vxc_sp1()

# Overwrite time coordinate to use fs instead of au as unit
au2fs = 0.024189
nx = nx.assign_coords(t=nx.t * au2fs)
nx.t.attrs["units"] = "fs"
vxc = vxc.assign_coords(t=vxc.t * au2fs)
vxc.t.attrs["units"] = "fs"

hv_nx = hv.Dataset(nx).to(hv.Curve, "x", "density_sp1").opts(ylabel="n(x)")
hv_vxc = hv.Dataset(vxc).to(hv.Curve, "x", "vxc_sp1").opts(ylabel=r"$v_{xc}(x)$")
layout = (hv_nx + hv_vxc).cols(1)
layout = layout.opts(title="e-H scattering computed at the ALDA level. | {dimensions}", sublabel_format = "", fig_size=250)
hv.output(layout, holomap="gif", fps=6)

Tutorial Validation Checks¶

import numpy as np
from tutorial_helpers.extract_peak import extract_peak

np.testing.assert_allclose(run.scf.info.get_total_energy(), -0.64659389, rtol=0.001)

x_coords = nx['x'].to_numpy()

np.testing.assert_allclose(extract_peak(x_coords, nx.sel(step=0), order=1),
                           (-10.05, 2.3223973528598822, 0.386521984833474), rtol=0.001)

np.testing.assert_allclose(extract_peak(x_coords, nx.sel(step=0), order=2),
                           (10.05, 3.725243372817074, 0.252187127109816), rtol=0.001)

np.testing.assert_allclose(extract_peak(x_coords, nx.sel(step=50), order=2),
                           (7.8, 4.209692330906179, 0.228832243413781), rtol=0.001)

np.testing.assert_allclose(extract_peak(x_coords, nx.sel(step=200), order=2),
                           (1.05, 7.009000523934921, 0.138043100711761), rtol=0.001)

np.testing.assert_allclose(extract_peak(x_coords, nx.sel(step=300), order=2),
                           (-3.45, 4.454593949246828, 0.100725782895539), rtol=0.001)