Absorbing boundaries¶
In this tutorial, we explore how to work with absorbing boundaries. For this, we will model the dynamics of an electron wavepacket scattering in a one-dimensional box and see how we can employ absorbing boundaries to efficiently avoid reflection at the edge of the simulation box.
Ground state¶
We will start with an input file containing only the description of a 1D hydrogen atom, modeled by a soft Coulomb potential. The details of this physical system are not really relevant here, as we will replace the electronic state with a user-defined state for the time-dependent calculation.
!mkdir -p 8-absorbing-boundaries/ground-state
%%writefile 8-absorbing-boundaries/ground-state/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 8-absorbing-boundaries/ground-state/inp
!cd 8-absorbing-boundaries/ground-state && octopus
These input variables should already be familiar as they have been explained in previous tutorials.
No absorbing boundaries¶
We are now ready to perform a simulation of the time-evolution of a wave-packet, without absorbing boundaries. For this, we use the following input file:
!cp -rT 8-absorbing-boundaries/ground-state 8-absorbing-boundaries/no-absorbing-boundaries
%%writefile 8-absorbing-boundaries/no-absorbing-boundaries/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
%Occupations
1 | 0
0 | 0
%
a = 0.1
x0 = 10.0
p = 5.
%UserDefinedStates
1 | 1 | 1 | formula | "exp(-a*(x-x0)^2+i*p*(x-x0))"
%
%Output
density
potential
%
OutputFormat = axis_x
TDPropagationTime = 0.5 * femtosecond
TDTimeStep = 0.01
Overwriting 8-absorbing-boundaries/no-absorbing-boundaries/inp
!cd 8-absorbing-boundaries/no-absorbing-boundaries && octopus
In order to introduce the electron wavepacket to scatter on the hydrogen atom, we have done the following change:
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=-10\) (a.u.) is the initial position of the wavepacket, and \(p\) its velocity, taken to be \(p=5\) a.u., such that the wavepacket is sent directly through the hard wall of the simulation box.
Note that we used a separate directory for the calculation so we can compare the data later with the case with absorbing boundaries.
While Octopus provides several methods, we will use the complex absorbing potential (CAP) method. This is done by adding:
AbsorbingBoundaries = cap
ABCapHeight = -1
%ABShape
40 | 50 | "abs(x)"
%
The variable AbsorbingBoundaries selects the type of absorbing boundaries. The height of the complex potential is controlled by ABCapHeight, and ABShape allows specification of the shape. Here we are absorbing in the region \(40<|x|<50\). In practice, the height of the complex potential, as well as the width of the absorbing region, need to be selected with care, as the reflection error (i.e. how much of the wavepack is not absorbed) is determined by these parameters. More details can be found in the implementation paper.
!cp -rT 8-absorbing-boundaries/ground-state 8-absorbing-boundaries/absorbing-boundaries
%%writefile 8-absorbing-boundaries/absorbing-boundaries/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
%Occupations
1 | 0
0 | 0
%
a = 0.1
x0 = 10.0
p = 5.
%UserDefinedStates
1 | 1 | 1 | formula | "exp(-a*(x-x0)^2+i*p*(x-x0))"
%
%Output
density
potential
%
OutputFormat = axis_x
TDPropagationTime = 0.5 * femtosecond
TDTimeStep = 0.01
AbsorbingBoundaries = cap
ABCapHeight = -1
%ABShape
40 | 50 | "abs(x)"
%
Overwriting 8-absorbing-boundaries/absorbing-boundaries/inp
!cd 8-absorbing-boundaries/absorbing-boundaries && octopus
The result of the calculation is shown in the figure below. Without absorbing boundaries, we observe an artificial reflection of the electron wavepacket at the border of the simulation box, which would cause a problem if we were interested in the electron dynamics inside our simulation box. As expected, this is fixed by using absorbing boundaries, which efficiently damp the out-going wavepacket, avoiding the use of extremely large simulation boxes, which are numerically very expensive.
from postopus import Run
import holoviews as hv
# Fetch density from run without absorbing boundary
run_no_absorb = Run("8-absorbing-boundaries/no-absorbing-boundaries")
density_no_absorb = run_no_absorb.td.density_sp1()
# Fetch density from run with the absorbing boundary
run_with_absorb = Run("8-absorbing-boundaries/absorbing-boundaries")
density_with_absorb = run_with_absorb.td.density_sp1()
# Overwrite time coordinate to use fs instead of au as unit
au2fs = 0.024189
density_no_absorb = density_no_absorb.assign_coords(t = density_no_absorb.t * au2fs)
density_no_absorb.t.attrs["units"] = "fs"
density_with_absorb = density_with_absorb.assign_coords(t = density_with_absorb.t * au2fs)
density_with_absorb.t.attrs["units"] = "fs"
# Now use holoviews to create a plot for each time step
hv_density_no_absorb = hv.Dataset(density_no_absorb)
curve_no_absorb = hv_density_no_absorb.to(hv.Curve, "x", "density_sp1", label="No Abs.")
hv_density_with_absorb = hv.Dataset(density_with_absorb)
curve_with_absorb = hv_density_with_absorb.to(hv.Curve, "x", "density_sp1", label="CAP")
# You can use one of the following methods to display the time series:
# --- slider controlled plot ---
#
#hv.extension("bokeh", "matplotlib")
#(curve_no_absorb * curve_with_absorb).opts(
# title="Effect of the absorbing boundaries on an outgoing electron wavepacket.",
# width=600, height=600,
#)
#
# ------------------------------
# ---- animated plot ----
hv.extension("matplotlib")
layout = (curve_no_absorb * curve_with_absorb).opts(
title="Effect of the absorbing boundaries on an outgoing electron wavepacket. ({dimensions})", fig_size=200,
)
hv.output(layout, holomap='gif', fps=5)
# ------------------------
Tutorial Validation Checks¶
from tutorial_helpers.extract_peak import extract_peak
import numpy as np
np.testing.assert_allclose(extract_peak(density_no_absorb['x'].to_numpy(), density_no_absorb.sel(step=0).to_numpy(), order=1),
(10.05, 3.7252435569845748, 0.252187127109816), rtol=0.001)
np.testing.assert_allclose(extract_peak(density_with_absorb['x'].to_numpy(), density_with_absorb.sel(step=0).to_numpy(), order=1),
(10.05, 3.7252435569845748, 0.252187127109816), rtol=0.001)
np.testing.assert_allclose(extract_peak(density_no_absorb['x'].to_numpy(), density_no_absorb.sel(step=500).to_numpy(), order=1),
(34.95, 6.39791692973688, 0.153770883463528), rtol=0.001)
np.testing.assert_allclose(extract_peak(density_with_absorb['x'].to_numpy(), density_with_absorb.sel(step=500).to_numpy(), order=1),
(34.95, 6.397915604360822, 0.153770871402382), rtol=0.001)
np.testing.assert_allclose(extract_peak(density_no_absorb['x'].to_numpy(), density_no_absorb.sel(step=1000).to_numpy(), order=1),
(40.2, 10.480840685545473, 0.0925532093404977), rtol=0.001)
np.testing.assert_allclose(extract_peak(density_with_absorb['x'].to_numpy(), density_with_absorb.sel(step=1000).to_numpy(), order=1),
(38.85, 10.014643418747326, 0.00182641037689572), rtol=0.001)
np.testing.assert_allclose(extract_peak(density_no_absorb['x'].to_numpy(), density_no_absorb.sel(step=1500).to_numpy(), order=1),
(15.75, 16.09403845154931, 0.0611844548534031), rtol=0.001)
np.testing.assert_allclose(extract_peak(density_with_absorb['x'].to_numpy(), density_with_absorb.sel(step=1500).to_numpy(), order=1),
(13.95, 13.190518933691418, 0.00126262769255595), rtol=0.001)
np.testing.assert_allclose(extract_peak(density_no_absorb['x'].to_numpy(), density_no_absorb.sel(step=2000).to_numpy(), order=1),
(-7.35, 22.235482142568546, 0.0439385570304855), rtol=0.001)
np.testing.assert_allclose(extract_peak(density_with_absorb['x'].to_numpy(), density_with_absorb.sel(step=2000).to_numpy(), order=1),
(-11.55, 16.57288220710913, 0.000990951583529114), rtol=0.001)